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ABSTRACT 


STABILITY OF THREE-DIMENSIONAL COMPRESSIBLE 

BOUNDARY LAYERS. 

Samarasingham Jeyasingham 
Old Dominion University, 1999 
Director: Dr. P. Balakumar 

A program is developed to investigate the linear stability of three-dimensional 
compressible boundary layer flows over bodies of revolutions. The problem is for- 
mulated as a 2D eigenvalue problem incorporating the meanflow variations in the 
normal and azimuthal directions. Normal mode solutions are sought in the whole 
plane rather than in a line normal to the wall as is done in the classical ID stability 
theory. The stability characteristics of a supersonic boundary layer over a sharp 
cone with 5° half-angle at 2° angle of attack is investigated. The ID eigenvalue 
computations showed that the most amplified disturbances occur around x 2 = 90° 
and the azimuthal mode number for the most amplified disturbances range between 
m ~ —30 to -40. The frequencies of the most amplified waves are smaller in the 
middle region where the crossflow dominates the instability than the most amplified 
frequencies near the windward and leeward planes. The 2D eigenvalue computations 
showed that due to the variations in the azimuthal direction, the eigenmodes are 
clustered into isolated confined regions. For some eigenvalues, the eigenfunctions are 
clustered in two regions. Due to the nonparallel effect in the azimuthal direction, 
the eigenmodes are clustered into isolated confined regions. For some eigenvalues, 
the eigenfunctions are clustered in two regions. Due to the nonparallel effect in the 
azimuthal direction, the most amplified disturbances are shifted to 120° compared 
to 90° for the parallel theory. It is also observed that the nonparallel amplification 
rates are smaller than that is obtained from the parallel theory. 

0 The format of this thesis is based on American Institute of Aeronautics and Astronautics Journal 
and was typeset in IAT]?X2 e by the author. 
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CHAPTER 1 
INTRODUCTION 


The transition of viscous flows from laminar to turbulent is one of the most 
challenging problems in fluid mechanics. In a typical commercial aircraft, 50 percent 
of the total drag force is due to the skin friction drag (Hefner and Bushnell [35] 1979). 
In supersonic and hypersonic vehicles, the equilibrium temperature on the surface 
determines the quality of material to be used in making the surfaces. Prediction of 
the transition onset point, the location at which the laminar boundary layer starts 
to become turbulent, is critical in the design of aerodynamic vehicles. The accurate 
prediction of the transition onset point is also very important in the application of the 
laminar flow control (LFC) methods to subsonic and supersonic aircrafts. Successful 
prediction of transition onset point depends on understanding the transition process 
and transferring the understandings onto a prediction tool. 

Though there are several mechanisms and routes to go from a laminar to 
a turbulent state, in quiet environments all of them, in general, follow these funda- 
mental processes. 

- Receptivity 

- Linear instability 

- Nonlinear stability and saturation 

- Secondary instability 

- Breakdown to turbulence. 

In the receptivity process, the unsteadiness in the environment and the in- 
homogeneities in the geometry generate instability waves inside the flow. In quiet 



environments , the initial amplitudes of these instability waves are small compared 
to any characteristic velocity and length scales in the flow. Goldstein [36] (1983 
a) theoretically explained using asymptotic methods how the Tollmien-Schlichting 
waves (T-S waves) are generated near a leading edge of a flat plate by the long wave- 
length acoustic disturbances and in a companion paper [37] (1985 b) described the 
scattering of T-S waves from the acoustic disturbances by the streamwise variations 
in surface geometries. In the second stage, the amplitudes of these instability waves 
grow exponentially downstream and this process is governed by the linearized Navier- 
Stokes equations. Further downstream, the amplitude of the disturbances become 
large and the nonlinear effects inhibit the exponential growth and the amplitudes 
of the disturbances eventually saturate or attain singular values. In the next stage, 
these finite amplitude saturated disturbances become unstable to two and/or three 
dimensional disturbances. This is called secondary instability and can be analyzed 
using Floquet theory, Herbert (1988) [38]. Beyond this stage the spectrum broadens 
due to complex interactions and further instabilities and the flow becomes turbulent 
in a short distance downstream. In this thesis, we investigated the linear instability 
of three dimensional supersonic boundary layers over a sharp cone at an angle of 
attack. 

Depending on the boundary layer profiles and flow parameters different 
types of linear instability waves are generated inside the boundary layers. For sta- 
bility analyses, boundary layers can be divided into following groups. 

Incompressible Flows Compressible Flows 

{ ^ f \ . 

Two-dimensional Three-dimensional Two-dimensional Three-dimensional 


In incompressible two-dimensional boundary layers, the disturbances may 
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Fig. 1.1 Schematic of velocity components within a swept-wing boundary layer illus- 
trating the definition of the cross-flow velocity (Reed & Sane). 

take one of the following instability waves. 

(1) inviscid instability (Rayleigh instability) 

(2) viscous instability (Tollmien-Schlichting instability ). 

Inviscid instability waves are generated by the interaction of inertia and 
pressure forces. This instability exists when the boundary layer contains an inflection 
point in the velocity distributions. The examples are boundary layer profiles in 
adverse pressure gradients, wakes, jets and in separated flows. Viscous instability 
waves are generated by the interaction of inertia, pressure and viscous forces. In this 
type of instability, energy is transferred from mean flow to the disturbance motion 
through the action of viscosity. Hence, these flows are inviscidly linearly stable but 
unstable to infinitesimal disturbances at finite Reynolds numbers. This instability 
is found in the plane Poiseuille flow and in boundary layer regions where pressure 
gradients are small, e.g, Blasius boundary layer. 
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Three-dimensional boundary layers are defined as the flows where the in- 
viscid streamlines are curved. In a three-dimensional boundary layer, there exists 
another instability mode, called cross-flow instability. When the inviscid streamlines 
are curved, there exists a pressure gradient in the direction normal to the inviscid 
streamlines. Inside the boundary layer, due to the viscous effects, the velocity is 
smaller than that in the inviscid region. Hence, this pressure gradient causes a veloc- 
ity component inside the boundary layer that is perpendicular to the inviscid- velocity 
vector. This component is called cross-flow. A schematic diagram showing differ- 
ent components of the velocity inside a boundary layer is given in figure 1.1. The 
cross-flow velocity profile has a maximum velocity somewhere in the middle of the 
boundary layer and goes to zero on the body surface and at the boundary layer edge, 
therefore exhibiting an inflection point. The description of the instability caused by 
the cross-flow was first given in the classic paper by Gergory, Walker, Stuart [40] 
(1956). When the cross-flow component is combined with the velocity component 
in the inviscid direction, they form a mean velocity profile that has an inflection 
point at which the velocity is zero. This permits a neutral disturbance with zero 
frequency. This neutral disturbance appears as vortices that all rotate in the same 
direction and take on the form of the familiar “cat’s eye” structure when viewed in 
the stream direction. This phenomena is observed in several flow geometries such 
as rotating cones ( figure 1.2 ), swept wings (figure 1.3), spheres (figure 1.4) and 
rotating disks ( figure 1.5). 

The stability characteristics of compressible two-dimensional and axi-symmetric 
boundary layers have been thoroughly investigated (Lees and Lin (1946) [41] , 

Lees and Reshotko(1962) [42], Mack(1969) [43], Gaponov(1981) [44], Malik and 
Spall(1991) [45]). The main conclusions that are drawn from these investigations 
can be summarized as follows. (1) The quantity (p'U)' = where p and U 

are the density and the streamwise velocity distributions of the meanflow and y the 
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coordinate normal to the wall, plays the same role in the compressible stability the- 
ory as U" does in the incompressible theory. The location where (p'U)' = 0 is called 
the "Generalized inflection point”. In most of the compressible boundary layers, the 
density profile has an inflection point p" = 0 and due to this compressible boundary 
layers exhibit generalized inflection points. As a consequence, the flat-plate/ cone 
compressible boundary layers are unstable to purely inviscid disturbances. This is 
one of the important difference between the instability of incompressible and com- 
pressible flows. (2) When the mean flow relative to the neutral disturbance phase 
velocity becomes supersonic over some portion of the boundary layer, there exist sev- 
eral unstable modes. For two-dimensional disturbances, it is the first of the additional 
modes, called the second mode, which is the most unstable at all Mach numbers for 
which the relative flow is supersonic. (3) At higher Mach numbers, the viscosity has 
only a stabilizing influence on the boundary layer. (4) Considering three-dimensional 
disturbances, the amplification rate of the first mode increases while the amplification 
rate of the higher modes decreases with the increasing waveangles. In the inviscid 

limit, the phase velocity of the first mode varies from 1 - h to where M ,S the 
free stream Mach number and is the mean flow velocity at which the generalized 
inflection {p'U)' = 0 occurs. To have unstable first mode disturbances the general- 
ized inflection point should appear above the mean velocity 1 " Jg- For the second 
mode, the phase velocity C varies from C, to 1. At low Mach numbers there exists 
no supersonic region near the wall relative to the phase velocity C„ hence no second 
mode instability exists at low Mach numbers. The supersonic region first starts to 
appear in the inviscid limit at M = 2.2 in an insulated flat plate boundary layer. 
The lowest Mach number at which the unstable second mode region has been located 
at finite Reynolds numbers is M = 3. The second mode instability increases with 

increasing relative supersonic region. 

Linear stability of axi-symmetric three-dimensional compressible boundary 
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layers were investigated by Balakumar and Reed (1991) [19]. As discussed pre- 

viously, at low Mach numbers, in compressible flows the most amplified waves are 
oblique while in incompressible flows they are two-dimensional. This is due to the 
fact that in supersonic flows, the amplification rate increases with decreasing Mach 
numbers. Since in an oblique direction, the effective Mach number decreases hence 
the amplification rate increases for three-dimensional waves. For a free-stream Mach 
number of M = 3, the wave angle of the most amplified wave is inclined at about 
55°. As described earlier, in incompressible flows the cross- flow velocity component 
introduces a new instability called “cross-flow instability” and the disturbances are 
inclined very close to the cross-flow direction. In compressible flows, the cross-flow 
basically increases the amplification rates of the first mode and makes the most am- 
plified disturbances inclined more towards the cross-flow direction. Balakumar and 
Reed’s calculation showed that the amplification rate of the first mode is increased 
by a factor of 2 to 4 due to the cross-flow compared with a two-dimensional flow and 
this increase decreases with increasing Mach number. It was also shown that the 
waveangles of the most amplified waves are increased by about 10° and the effect of 
the cross-flow on the second mode is as expected small. 

In general three-dimensional boundary layers, the mean boundary layer 
profiles vary in all three directions: streamwise , spanwise or azimuthal and normal 
directions. However, in the high Reynolds number boundary layer flows the varia- 
tions in the streamwise and in the spanwise directions are smaller than that in the 
normal directions. In the classical stability theory, these variations in the streamwise 
and in the spanwise directions are neglected and it is assumed that the boundary 
layer profiles vary only in the normal direction and uniform in the other two direc- 
tions. This makes the coefficients of the linearized Navier-Stokes equations to be 
independent of the streamwise and spanwise coordinates and permit one to seek a 
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solution in the normal mode form 




(u: 


where x, , ,r 2 and * 3 are the coordinates in the streamwise, spanwise and in the normal 
directions, t is the time. $(x 3 ) is the eigenfunction , a,/? are the wavenumbers in 

directions and w is the frequency. In general, a,/? and 


the streamwise and spanwise 
u> are complex. 


a = + ioi , 

/3 = f3 r + ifti ? 

a? = u; r + * 


(1.2) 

(1.3) 

(1.4) 


When the normal mode form Eq. 1.1 is substituted into the linearized Navier-Stokes 
equations , a homogeneous system of ordinary differential equations with homoge- 
neous boundary conditions are obtained. The solution of which yields a dispersion 
relation among a, /? and ui of the form 


F{aJ,u>) = 0- 


( 1 . 5 ) 


The real and imaginary parts of the relation yield two equations in terms 
of the six unknown parameters (a„a,) , (A. ft) :ultl (“r.w.) • To determine these 
six unknowns one needs to specify four additional conditions. Two approaches are 
generally used to overcome these difficulties. One is called the temporal eigenvalue 
approach in which the wavenumbers o,/3 are prescribed as real numbers and the 
complex frequency u is solved from the dispersion relation Eq. 1.5. The other method 
is called spatial eigenvalue approach in which the frequency tv and the spanwise 
wavenumber 0 are prescribed as real numbers and the complex wavenumber a is 

solved from the dispersion relation. 


8 


Temporal Problem 

a ~ a r , real prescribed 

0 = d r , real prescribed 

u! = lj t + iijJi , solved from eigen relation 
uj r = frequency of the disturbance 
uji = amplification rate of the disturbance in time 
cu'i < 0 , the boundary layer is stable 

= 0 , the boundary layer is neutrally stable 

> 0 , the boundary layer is unstable 


Spatial Problem 

uj = uJ r , real prescribed 

0 = 0 T y real prescribed 

a = a r + ic*i , solved from eigen relation 
a r — wavenumber in the streamwise direction 
—a,* = growth rate in the streamwise direction 
—a, < 0 , the boundary layer is stable 

= 0 , the boundary layer is neutrally stable 

> 0 , the boundary layer is unstable 

In the temporal method, lJ{ measures the amplification of the disturbances 

in time and in the spatial method, -a, measures the growth rate of the disturbances in 
the streamwise direction. The boundary layer is stable, neutrally stable or unstable to 
small disturbances depending on whether the amplification rateu;, or — a x is less than 
, equal to or greater than zero respectively. Most of the linear stability computations 
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have been performed based on these classical approaches. In this work, this approach 
is called ID method. 

The next step is to estimate the correction to the eigenvalues that is ob- 
tained from the parallel theory due to the small variation of the meanflow in the 
streamwise and in the spanwise directions. Three methods are available to compute 
the evolution of small disturbances in a non-parallel flow. One and the oldest method 
is the multiple scale approach (Saric , Nayfeh [46]). The second is the Parabolized 
Stability Equations (PSE) approach (Herbert 1979 [53] )and the third method is 
to solve the full Navier-Stokes equations in a non-parallel flow ( Fasel [49] . Joslin 
[51]). In the linear regime, instead of solving the full Navier-Stokes equation, the 
linearized Navier-Stokes equations are solved [52]. In the multiple scale and in the 
PSE methods the disturbances are written in the form 


$(x 1 ,x 3 ) = $ 0 (x 1 ,x 3 )e ! -/'' l(xi)[ixi+, ^ 2 - ,w< • (1.6) 

Here u; is the real frequency, ft is the real spanwise wavenumber , a(xi ) is the stream- 
wise wavenumber which is a function of x\ and $ 0 (xi,X 3 ) is the amplitude function 
which is a function of both X\ and x$. This form of the representation is mathemat- 
ically and physically meaningful in a meanflow which varies only in the streamwise 
(xi) and normal (£ 3 ) direction and is uniform in the spanwise (x 2 ) direction. These 
approaches are used to compute the evolution of the disturbances in two-dimensional 
Blasius type boundary layers and in quasi-three dimensional, infinite swept wing, 
boundary layers ( El Hady [48], Herbert [33] and Malik [34] ). 

The objective in this work is to investigate the stability and the evolu- 
tion of disturbances in fully three-dimensional boundary layers. By the fully three- 
dimensional boundary layers it is meant boundary layer flows over finite wings, 
flow over non-axisymmetric geometries like ellipsoids, delta wings and flow over axi- 
symmetric geometries at angles of attack etc. Specifically, in this thesis the stability 
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of supersonic boundary layer flows over a sharp axi-symmetric cone at an angle of 
attack (King 1992 [54]) is investigated. In this case, the meanflow varies in all three 
directions, streamwise, azimuthal and in normal directions. Since the mean flow 
varies in the azimuthal direction, it is not possible to decompose the disturbances 
as a sum of Fourier components as is done in the axi-symmetric or in the infinite 
swept wing flows. The expectation is that since the instability is directly related to 
the local mean flow conditions, the eigenfunctions will be confined to a region in the 
azimuthal direction. Hence, the normal mode is written in the form 

ar 2 , a^ 3 ) = $ 0 (x, <x 2 ,x 3 )e‘f a ^ d *'- twt . (1.7) 

The azimuthal variation is included in the amplitude part $(x 1 ,X 2 ,X 3 ) 
which now becomes a strong function of x 2 and x 3 and a slowly varying function 
of xi. If the xi dependence of $ 0 and a are dropped, one can obtain an eigenvalue 
problem for a or w and 4 > 0 (:e 2 , £ 3 ) which is a function of x 2 and x 3 coordinates. This 
is called as 2D eigenvalue problem. 

Several experiments were performed to understand the stability and transi- 
tion of supersonic and hypersonic two dimensional and three dimensional boundary 
layers. The experiments can be divided into two groups. One is transition experi- 
ments in which the transition onset is measured at different flow conditions (Potter 
1974 [18], Krogmann 1977 [55], Stetson 1981 [9], King 1992 [54]). In these exper- 
iments, the unit Reynolds number effects and the effects of the angles of attack on 
the transition front are investigated. 

Since there is no length scale in flows over sharp cones, the transition 
Reynolds number should not change with the free stream unit Reynolds numbers. 
However, experiments performed in different wind tunnels at different Reynolds num- 
bers change with the free stream unit Reynolds number. Though transition is influ- 
enced by several factors, bluntness, angle of attack, vibration of the model, roughness, 
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free stream disturbance levels, spectral content of the free stream disturbance, the 
unit Reynolds number effects are generally attributed to the radiated acoustic field 
from the boundary layers on the tunnel walls. Main conclusion is that in quiet tun 
nels at low unit Reynolds numbers the boundary layers on the wall remain laminar 
and the effect of unit Reynolds number on the transition Reynolds number is mini- 
mal. In noisy wind tunnels the boundary layers on the walls become turbulent and 
the transition Reynolds numbers decrease with the unit Reynolds numbers. 

The experiments performed over sharp and blunt cones at small angles of 
attack show that the transition front moves downstream on the windward ray and 


moves downstream on the leeward ray. 


The second type of experiments are the stability experiments where the 
stability characteristics of the boundary layers are investigated. Stetson and his co- 
workers systematically investigated the stability characteristics of hypersonic bound- 
ary layers in natural conditions (Stetson et.al., 1983 [10], 1984 [11], 1985 [M> 1986 
[13], 1989 [14]). Their results are summarized in a review paper by Stetson and 

Kimmel(1992) [15]. Experiments were performed on sharp and blunt cones with i° 
half angle at a free stream Mach number of 8 at zero and nonzero angles of attack 
and with adiabatic and cooled surface conditions. 


The experiments agree with the theoretical predictions that the transition 
in hypersonic boundary layers is dominated by the high frequency second mode type 
disturbances. However, the measured growth rates are much smaller than that is 
computed from the linear stability theory. It is also observed that the small nose 
bluntness increases the critical Reynolds number from that for a sharp cone. Flow 
over a sharp cone at an angle of attack showed that the boundary layer along the 
windward direction becomes more stable and on the leeward side it becomes more 
unstable compared to that is obtained at zero angle of attack. Hence transition on 
the leeward side moves downstream and in the leeward side moves upstream from 
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that for a sharp cone. The frequencies of the most amplified disturbances on the 
windward side are larger than that is obtained for a sharp cone at angle of attack. 
This is attributed to the finding that the frequency of the most amplified waves in 
hypersonic boundary layers depends inversely on the boundary layer thickness which 
is smaller than that is obtained at zero angle of attack, hence the frequency of the 
most amplified waves are higher. Unit Reynolds number effects are investigated by 
measuring the spectral content of the free-stream disturbances and the range of the 
most amplified frequencies for the the boundary layer. It is observed that most of the 
energy in the freestream disturbances are contained in the low frequency disturbances 
and the spectrum for the high frequency disturbances are very small. But it was 
concluded that the high frequency disturbances with very small amplitudes are are 
sufficient to initiate the second mode disturbances. It was also observed that if the 
frequency of the most amplified wave for the boundary layer is much higher than 
the frequency limit in the free-stream disturbance spectrum, then the most amplified 
wave will not be initiated and the transition is dominated by the smaller frequency 
waves which are excited by the free-stream disturbances. Instantaneous structure and 
the ensemble-averaged structure of the second mode instability waves in a hypersonic 
boundary layer was measured by Kimmel et. al. (1997) [22], Poggie(1997) [7] 

in natural conditions. It is observed that the second mode disturbances travel as 
wavepackets confined to a small region in the circumferential direction. Recently 
Poggie et. al. (1998) [6] investigated the stability and transition of a hypersonic 

three-dimensional boundary layer over an elliptic cone. Transition front appears 
asymmetric with early transition close to the minor axis and delayed transition close 
to the major axis which is similar to that is observed in the flow over a sharp cone 
at an angle of attack. 

King (1992) [54] investigated the transition of a three-dimensional bound- 
ary layer in NASA Langley’s Mach 3.5 quiet tunnel. The experiments were performed 



for a 5° half angle sharp cone at various angles of attack. The transition is dominated 
by crossflow dominated instability and the transition onset point moves downstream 
near the windward side and moves upstream near the leeward side. This case is cho- 
sen to investigate the stability characteristics of a three-dimensional boundary layer 
using the 2D eigenvalue method. 



a) Rotational speed below critical- no instability apparent 




lly . 
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b) Rotational speed slightly above critical - spiral streaks are observed 



c) Rotational speed far above critical- spiral streaks are observed, then 
secondary instability superposed on the vortices, then transition 

Fig. 1.2 Flow visualization for a spinning cone (Kobayashi et al. 1983) 


Fig. 1.3 Naphthalene surface patterns illustrating cross-flow vortices (Saric & Yeates 
1985). Flow is from left to right. 



Fig. 1.4 Flow visualization for a rotating sphere: spiral streaks are observed , then 
the secondary instability superposed on the vortices, then transition (Kohama &; 
Kobayashi 1983b). 



Fig. 1.5 Flow visualization illustrating the spiral vortices on a disk (Kohama 1984a). 
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CHAPTER 2 
FORMULATION 


This chapter describes the formulation of the linear stability theory of three- 
dimensional compressible boundary layer and the numerical procedure that is used 
to solve these equations for a model case - supersonic flow over a sharp cone at an 
angle of attack. The first section deals with the derivation of the linear Parabolized 
Stability Equations in its most generic form i.e, in generalized curvilinear coordinate 
system. In the second section, the classical one-dimensional stability equations and 
the two-dimensional stability equations are derived from the general theory as special 
cases. 

As described in the introduction, the stability problem is formulated in two 
different methods. The first methodology is to make a locally parallel flow assump- 
tion i.e, neglect the meanflow variations in the streamwise and azimuthal directions 
and formulate the problem as a ID eigenvalue problem. In other words, the solutions 
to the stability equations are sought in a line normal to the cone surface. In fact, 
large portion of the literature on stability computations of three-dimensional and two 
dimensional boundary layers are performed as ID eigenvalue problems. But the par- 
allel flow assumption will not hold true for flows which are highly three-dimensional 
and the stability computations and transition prediction will not be correct. Hence 
in an attempt to construct a near approximate solution to the full Navier-Stokes 
equation, the variation of the meanflow in the azimuthal direction is incorporated 
in the second method and formulated as a 2D eigenvalue problem. Thereby normal 
mode solution is sought in a plane at a streamwise location, rather than in a line as 
in the former method. In both cases, ID and 2D eigenvalue approaches, the problem 
can be formulated as temporal or a spatial stability methods. In the temporal stabil- 
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itv the wave number in streamwise direction a is known and the desired eigenvalue 
is the temporal frequency a;, whereas in spatial stability u; is fixed and a is sought. 

There are two classes of numerical methods that can be used for the com- 
putations of the temporal or spatial eigenvalues from the eigenvalue problems: global 
and local methods. In the global method a generalized eigenvalue problem is set up 
and the eigenvalues are obtained using some standard algorithms such as QZ avail- 
able in the public-domain software library LAPACK. Here, a guess for the eigenvalues 
is not required. On the other hand, in the local method a guess for the eigenvalue 
is required and only the eigenvalue that happens to be in the neighborhood of the 
guessed value is computed using some iterative techniques. These methods will be 
discussed in detail in the last section of this chapter. 


2.1 Formulation of the Stability Theory 

The growth or decay of infinitesimal perturbations superimposed on a laminar flow 
is the focus of linear stability theory. The linear stability theory analyzes the char- 
acteristics of the instabilities of the mean laminar flow over the surface of interest. 
Transition prediction is basically composed of two tasks; (1) accurate calculation of 
the viscous flow field over the the body, (2) calculation of the amplification rate. The 
stability properties of two-dimensional incompressible and compressible boundary- 
layers and three-dimensional incompressible and compressible boundary-layers were 
discussed in the first chapter. In this section, the derivation of the linear stability 
equation, starting from the compressible three-dimensional Navier-Stokes equations 
in orthogonal curvilinear coordinate system is presented. Normal mode method is 
chosen in the formulation. Since the final expression for the linear stability equa- 
tions in the generalized curvilinear coordinate system contain numerous terms only 
linearized continuity equation will be given. 
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Fig. 2.1 Orthogonal curvilinear coordinate system 


2.1.1 Governing Equations 

The derivation of the stability equations starts from the Navier-Stokes equations for 
three-dimensional compressible flow for an ideal viscous gas in orthogonal curvilinear 
coordinate system . The governing equations are the continuity equation, the mo- 
mentum equations in the streamwise, azimuthal and normal directions and energy 
equation. In addition, perfect gas equation of state is used to 'close’ the system 
of fluid dynamic equations. A body-oriented coordinate system is used as shown 
in figure 2.1, with £i taken along the streamwise direction, X 2 along the azimuthal 
direction and X 3 normal to the surface. The Navier-Stokes equations in the vector 
form are: 

Continuity equation 


dp 

dt 


+ V-pV - 0 


( 2 . 1 ) 
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Momentum equations 

P^X + V. VVJ =pf + v n tj ; i,j = 1,2,3 ■ (2.2) 

Energy equation 

P c v + v • vr ) = v • ( KVr) + $ ' (2,3) 


Equation of state 


p = pRT • 


(2.4) 


The first term on the right hand side of the momentum equations, Eq. 2.2, 
is the body force per unit volume. Ilq represents the component of the stress tensor, 
which consists of normal stresses and shearing stresses. $ in the energy equation, 
Eq. 2.8, is the dissipation function. 

The expanded form of these expressions as given by Anderson et.al. [32] is 
adopted here. 

Continuity equation 


+ 


1 


dp 

dt ' h\h 2 h 3 


-^—{h 2 h 3 pui) + -^—{hzhipui) + -x—{hih 2 pu 3 ) 
dxi ox 2 ox 3 


= 0 


(2.5) 


aq momentum 
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dt 
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h 2 dx 2 
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_d_ 
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dx 3 


1 dhi , n 1 dhi n 
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hih 2 dx 2 

1 dh 3 
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1 dh 2 
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x 2 momentum 



pu\ du 2 
hi dx x 


pu 2 du 2 
h 2 dx 2 


PU3 du 2 
h 3 dx 3 


pu\ dhi 
h\h 2 dx 2 


+ 


pui u 2 dhi 
hih 2 dx 2 


pu^u^ dh 2 pu\ dh 3 

h\h 3 dx 3 h 2 h 3 dx 2 


1 

h\h 2 h 3 


{h 2 h 3 n XiXi ) + — — (/jiA 3 n 


dxi 


dx 2 


% 2 X 2 


) "E {h\ h 2 Yl X 2 x 3 ) 


+n 


12^3 


1 dh 2 

h 2 h 3 dx 3 


n 


Xl x 2 


1 dh 2 

h \ h 2 c)x i 


^ 2 : 3 x 3 


1 dhz 

h 2 h 3 dx 2 


-n 


XlX! 


1 dhi 

hih 2 dx 2 


(2.6) 


a : 3 momentum 



pu i du 3 
hi dxi 


+ 


pu^du^ 
h 2 dx 2 


PU3 duz 
h 3 dx 3 


pu\ dhi pu\ dh 2 

hih 3 dx 3 h 2 h 3 dx 3 


puiu 3 dhz 
hih 3 dxi 
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^ 02^3 Obi 
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1 
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d 
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1 dh 3 
h\ h 3 dxi 
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1 dh 3 
h 2 h 3 dx 2 


n XlJ:i 


1 dh x 
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X 2 X 2 


1 dh 2 
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Energy equation 


pC v 


dT dT uidT_ u 2 dT u 3 dT' 
dt P dt hi dxi h 2 dx 2 ^ h 3 dx 3 


(2.T) 
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+ 7 — f~ L 7T — (h 2 h 3 Ui) + — — [hih 3 u 2 ) + 77 (hih 2 u 3 ) 

h \ /i 2 ^3 C/ 3 Tl < 7 X 2 C/X3 


1 ~ f) $ 

~j 7 7 77 ( Ar/i-2 ^ 3 ^x 113 ) + 77 ( /c h \ h 3 Ylx2x^ , 

h\h 2 h 3 ox\ ox 2 


+77 — (A-'/ii^n^j) + $ • (2.8) 

OX 3 

Here u\, u 2 and u 3 are velocities along the streamwise(.T!) , azimuthal(x 2 ) and 
normal(:r 3 ) directions respectively and p is the density, p is the pressure, p is the 
coefficient of viscosity and k is the coefficient of thermal conductivity, hi, h 2 and h 3 
are the metric coefficients along the coordinates x 2 , x 2 and x 3 respectively. In the 
generalized curvilinear coordinates, the dissipation function becomes 

$ = p 2(en 1 2 + e 22 2 + e 33 2 ) + e 23 2 + ei 3 2 + ei 2 2 — ^( e n + e 22 + e 33) ■ (2.9) 

where the expressions for the strains are 


&U — 


1 dui u 3 dh\ u 2 dhi 
hi dx x h x h 3 dx 3 hih 2 dx 2 


1 du 2 u 3 dh 2 u\ dh 2 
622 h 2 dx 2 **" h 2 h 3 dx 3 + hih 2 dxi 


1 du 3 ui dh 3 u 2 dh 3 
633 h 3 dx 3 hih 3 dxi + h 2 h 3 dx 2 


1 du 3 1 du 2 u 2 dh 2 u 3 dh 3 
623 h 2 dx 2 h 3 dx 3 h 2 h 3 dx 3 h 2 h 3 dx 2 


1 du 3 1 du 3 ui dhi u 3 dh 3 

613 h 3 dx\ hi dxi hih 3 dx 3 hih 3 dxi’ 
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1 du 2 1 dui 
e 12 — T" — h 


n 2 dh 2 ui dhi 


h\dx i h 2 0x 2 hih 2 dxi hih 2 dx 2 

The component of stress tensors appearing in the momentum 

Eq. 2.6 and Eq. 2.7 can be expressed as 

Pj _ 2 2-du.i 2 u 2 dhi 2 a 3 dh x 1 du 2 

' riXl P 3^ hi c>zi h x h 2 dx 2 + hih 3 dx 3 h 2 dx 2 


equations Eq. 2.6, 

^3 Oh 2 
h 2 h 3 dx 3 


u\ dh 2 1 du 3 U\ dh 3 u 2 dh 3 
hih 2 dxi h 3 0 x 3 hih 3 0 xi h 2 h 3 dx\ 


E-X1X2 — h 


1 du 2 
hi dxi 


xi 2 dh 2 
h x h 2 8 xi 


+ 


i du i 
h 2 0 x 2 


u x dh x 
hih 2 dx 2 ' 


n^ixa — P 


1 dui 
h 3 dx 3 


ui dh x 
h x h 3 dx 3 


+ 


1 du 3 

hi dx j 


U 3 dki 
hih 3 dx i 


2 

11x2X2 = ~P + r/* 


2 du 2 2 u 3 dh 2 2«i dh 2 

I ; -\ i 

h 2 dx 2 h 2 h 3 dx 3 h x h 2 dxi 


1 dui 
hi dxi 


u 2 dh\ 
hi h 2 dx 2 


u 3 dh x 1 du 3 iti dh 3 u 2 dh 3 
h x h 3 dx 3 h 3 dx 3 h x h 3 dxi h 2 h 3 dx 2 


( 2 . 11 ) 


n 


X2X3 


ii 1 du 3 
Re h 2 dx 2 


u 3 dhs 1 du 2 u 2 dh 2 
h 2 h 3 dx 2 h 3 dx 3 h 2 h 3 dx 3 \ 


n 


^3^3 


2 

= - p+ r 


2 du, 
h 3 dx 3 + 


2u\ dh 3 2u 2 dh 3 

h x h 3 dx 1 h 2 h 3 dx 2 


1 dui 
hi dxi 


u 2 dh\ 
h\ h 2 dx 2 


U3 dhx 1 du 2 u 3 dh 2 u\ dh 2 
h x h 3 dx 3 h 2 dx 2 h 2 h 3 dx 3 hxh 2 dxx 

2.1.2 Linearization of the equations 

The principle of classical stability theory evolves around the concept of determining 
whether a small disturbance introduced into a laminar boundary layer will amplify 
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or decay. If the disturbance decays, the boundary layer is stable and if the distur- 
bance grows it is called linearly unstable. In stability theory, the first step in the 
methodology of analyzing the evolution of small disturbances is to assume the * total 
flow' as composed of mean quantities and small disturbance quantities. 


% 2 •> £3 •> i (2.12) 

and in components form 

U 1 — f'i + wi 1 

l 1 2 = U 2 + ^2 1 
U 3 = l>3 + 11 3 , 

P = P + 7T , 

T = T + 6 , 

p = p + p , (2.13) 

where Q is the total quantity, Q is the mean quantity and q is the disturbance 
quantity. 

The stability equations are derived as follows: first, the expressions for 
total flow quantities, Eqs. 2.13, are substituted into the Navier-Stokes equations 
Eqs.( 2.6 - 2.8). Since it is assumed that the mean-flow terms satisfy the steady 
Navier-Stokes equations, the mean terms can be subtracted out, resulting in terms 
consisting of products of mean-flow and disturbance quantities (Q q) and the products 
of disturbance terms ( qq ). Secondly, in the linear theory, since the nonlinear terms, 
the products of infinitesimal disturbances, are of lower order than rest of the terms, 
they are neglected. Substituting the normal mode form for the disturbances, Linear 
Stability Equations (LSE) are obtained. 
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2.1.3 Nondimensionalization of the Equations 


As is customary and convenient, the LSE are written in nondimensional form using 
some characteristic quantities. The characteristic velocity is U 0 , which is taken to be 
the value of streamwise velocity at the edge of the boundary- layer in ID formulation 
and the boundary layer edge velocity at a reference station x 2 = 90° in the case of 
2D formulation; the characteristic length is L, which is given by the expression 


L = 



(2.14) 


The thermodynamic quantities are nondimensionalized by their corresponding boundary- 
edge values i.e, the characteristic density, temperature and molecular viscosity are 
p e , T e and p e respectively, the characteristic pressure is p e U 0 2 and the characteristic 
time is jj-. 

The non-dimensional quantities are 


«i 


7T = 


U 1 

Uo" 


IT 


Pe* Uo* 

TT~ W 

U 0 * 


(-'I T T + 1 


_ T 
T = — 
TS’ 


X\ = 


Xj_ 

L ’ 


u 2 


U 2 

u: ’ 

0 = fL 

T*' 


U 2 = rr „ 


X 2 = 


lh_ 

U 0 - 

T_ 

Pe*' 

x j1 

L ’ 




u 3 

U 0 * 


7T = K 
3 Do" 

_ T 

P = — , 

Pe 


X 3 


X3 

L ’ 


(2.15) 


Here the superscript * denotes the dimensional quantities. Some of the non-dimensional 
parameters which appear in the equations are defined below. 

Free-stream Mach number 


Up 


M = 


(2.16) 
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Prandtl number 

Pr = C p - ■ (2.17) 

AC 

Here Prandtl number is assumed to be constant and taken to be 0.7 all through the 
computations. The constants of specific heat at constant pressure and volume are 
related by 

CV = 7 CV (2.18) 

The values of ratio of specific heats 7 and gas constant R are taken to be 


7 = 1.4, (2.19) 

R = 287m 2 fs 2 °K- (2.20) 

2.1.4 Reduction of the Number of Variables 

The Linear Stability Equations contain perturbation terms, Ui, 112,113, 7r,0,p,fi and 
ac, which are the unknown variables. Since some of these variables are related to 
other variables by simple equations such as Eq. 2.4, Eq. 2.17 one can easily elimi- 
nate them from the stability equations. This will result in substantial reduction in 
computational effort and storage requirements. 

The coefficient of viscosity ft is assumed to be a function of temperature 
only and Sutherland’s viscosity law is used in the computations. 



( 2 . 21 ) 


where 


C = U0AK/T e , 


( 2 . 22 ) 


and T e is the temperature at the edge of the boundary-layer. Hence, the total vis- 
cosity term which is the sum of mean and perturbation quantities, can be written 



27 


as 


P + /< = f*{T + 0 ), 

and using Taylor series expansion for small temperature fluctuation one can write 


mt + f>) = /.m + + 


- „ 


Hence the disturbance p can be written as 


(2.23) 


dp Q 

u = — =y , 

^ dr 

using equation Eq. 2.24 the derivatives of p become 


(2.24) 


dpd6_ 

dxi fjT 2 dxi + dT dxi 

For reasons that would become clear later in this chapter it is chosen to 
express it in terms of density fluctuation p and temperature fluctuation 6 , rather 
than eliminating p from LSE. The equation of state in non-dimensional form is 


PjM 2 = pT , 


(2.25) 


The pressure fluctuation is related to density and temperature fluctuations by the 
expression 


7T 7 M 2 = pT + pd , 


(2.26) 


or. 


7 r = 


( 7 M 2 ) P+ ( 7 m 2 1 9 ‘ 


And the derivatives of t with respect to the coordinates are 


(2.27) 
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dir 
Ox i 


7A/ 2 / dxi \ 7 A/ 2 


dO 1 fdT\ 1 

dxi + tM 2 J p + 7M 2 



(2.28) 


Hence, the dependent variables in the stability equations are q = {iq, u 2 , u 3 . p, 0} T . 


2.1.5 Introduction of Harmonic Disturbances 


As discussed previously, the stability of a three-dimensional boundary layer is inves- 
tigated by seeking a solution of the form: 


q(xi,x 2 ,x 3 ,f) = ^(x^x^x^e 1 / 0 ^ lwt -f c . c . 


(2.29) 


Here q = (tq, u 2 , u 3 , p, 0} , is a complex amplitude function of the disturbance 
variables; a(x x ) is the wavenumber in the streamwise direction and uj is the temporal 
frequency of the disturbances. 

The first and second derivatives of q(aq,x 2 ,x 3 , t) with respect to the stream- 
wise coordinate X\ are 


dq 

dx 1 


iot(xi )q(xj, x 2 ,x 3 ) + 


dq 

dx 


S' 


J Otdx I - 


■tuft 


(2.30) 


d 2 q 
d^ 


-a 2 (x 1 )q(xi,x 2 ,x 3 ) + i ~q(xi, x 2 , x 3 ) 

UX\ 

+ 2ia~ + e'f aixi)dxi ~ iu,t ■ 

OX\ UX\ J 


(2.31) 


Since the variation of the amplitude part of the disturbance q in the stream- 
wise direction is very small, the term is neglected and the equation Eq. 2.31 
becomes 


d 2 q f 2 . da . dq 

~ (xi)q(xx,x 2 ,x 3 ) + *— q(xi,x 2 ,x 3 ) + 2 


j i fa(xi )dx 1 


-tu/t 


(2.32) 
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By substituting these expressions into the linear stability equations Eq. 2.6 - 
2.8 one can obtain a set of partial differential equations for q. These sets of equations 
incorporate the effect of meanflow variation in streamwise direction and approximate 
the full Navier-Stokes equations. The concept was first introduced by Herbert [33] 
and the resulting equations are termed Parabolized Stability Equations(PSE). 


The PSE for compressible flow in generalized coordinates contain number 
of terms, and hence, in order to save space, only the linearized continuity equation 
in non-dimensional form is given here. 

dp \ 


—tup + 


+ 


+ 


+ 


+ 


h\ h 3 dx\ 

1 dh x 
h x h 2 dx 2 

1 dh2 
h 2 h 3 dx 3 

1 dh 3 _ 


pUx + 


pU2 + 


pU 3 + 


h\ h 3 Ox i 
1 dhi 


pu\ + 


1 dh 2 _ 

+ -rr7r~P U3 + 


1 

dh x 

h\h 2 

dx 2 

1 

dh 3 

h 2 h 3 

dx 2 

1 

dhx 

h\h 3 

dx 3 

1 

dh 2 . 

h\h 2 

dx x 

1 

dh 3 

h 2 h 3 

dz 2 

1 


h\ /13 

5 ar 3 


plh + 


pU 2 + 


v x ^rx p ^K\ apJr dxx 


p+ 


h 2 dx 2 h 2 dx 2 


— , 1 dU 3 , U 3 dp 

pU 3 + T""X p + - 7 — -X 

h 3 dx 3 h 3 dx 3 

1 dp 


pu x + 


, o u i T 7 ~(* qu i + (2.33) 

h\ ux\ h\ ux\ 


l dp p du 3 

pU3 + T~T. U 3 + —J 

h 3 dx 3 h 3 dx 3 


= 0 ■ 


The sets of linear PSE equations for q can be written in matrix form as 


A 2 


d 2 q 

~0x 2 2 


+ ^3 


d 2 q 

d ^ 


+ ^-12 


d 2 q 
dx x dx 2 


+ ^13 


<9 2 q 

dx x dx 3 


+ ^23 


a 2 q 
dx 2 dx 3 
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+ Bl |i + Bi |i + ft |a + rq = o 

dxi OX 2 dx 3 


(2.34) 


Here, A 2 A 3 , Ai 2 , T 13 , A 23 , 5i, B 2 , B 3 and C are (5x5) complex matri- 
ces which are functions of mean-flow quantities and their derivatives, o, and the 
matrices h\, h 3 and their derivatives. 


2.1.6 Boundary Conditions 

At the solid wall, no-slip conditions apply to the disturbance velocities. The temper- 
ature perturbations are assumed to vanish at the solid boundary. This is a reasonable 
assumption since for almost any frequency of the gas, the temperature fluctuation 
will not penetrate into the solid boundary due to the thermal inertia of the solid 
body i.e, the wall can only remain at its mean temperature. 


txj = u 2 — u 3 = 0 = 0; at x 3 — 0 • (2.35) 

The boundary conditions in the far field are that disturbances decay to zero. 

«x, u 2 , u 3 , 6 — *• 0; as .r 3 — ► oo. (2.36) 

The numerical procedure which is used in the computations need one more 
boundary condition to be specified for the density perturbation. Since the den- 
sity does not vanish at the wall and the PSE equations are valid both at the solid 
boundary and at infinity, the linearized continuity equation or the normal momentum 
equation can be used as the fifth boundary condition at the wall and at the tarfield. 

2.1.7 Linear Stability Equations 


2D Probtejn 
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If the Xi derivative terms in the Parabolized Stability Equations are dropped, 
one can obtain 2D linear stability equations which are systems of partial differential 
equations in ar 2 and .r 3 . 


d 2 q 

W 


+ + A a — li- 

dx 3 2 dx 2 dx 3 


+ B. 


r)q 


+ f ?3 7 T~ h C q — 0 
ox 3 


(2.37 


This corresponds to applying classical parallel stability theory in the x\ 


direction and seeking normal mode solution of the form 


q{xi,x 2 ,x 3 ,t) = q(ar 2l x' 3 )e‘ Q(ri)i;i ,u " + c • c 


(2.38) 


Here a is the streamwise wavenumber, u> is the temporal frequency. In this work, 
the temporal stability computations are performed where 


a = a r , real prescribed , 

UJ — 0J r + * 


(2.39) 


(2.40) 


io is computed from the dispersion relation 

F(a,u}) = 0 • 


(2.41) 


/ D Pj^jblern 

In this method of formulation the variations of meanflow in the streamwise 
direction and azimuthal directions are neglected and the disturbance quantities are 
assumed to be periodic in the azimuthal direction. This is equivalent to applying 
classical stability theory in an and x 3 directions and to seeking normal mode solution 

of the form 

q(x u x 2 ,x 3 ,t) = q(x 3 )e'°*' +tm *>-'“ t +c-c - 


(2.42) 
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Applying the chain rule for the normal derivatives, one obtains 


dq(> 3 ) _ ( dr]\ dq 

dx 3 \dx 3 J dr]' 

d 2 q(x 3 ) ( drj \ 2 d 2 q / cl 2 r] \ dq 

dx 3 2 \dx 3 ) dr ] 2 \dx 3 2 J dr] 


(2.49) 


With these transformations the stability equations in the new coordinate // becomes 


A3 


dr] \ * d 2 q 
K dx 3 J drj 2 


+ 


im A 23 



+ B3 


dq 

dr] 


+ 


— m 2 A2 


d 2 77 


+ imB2 + C 


q = 0 


(2.50) 


The system of equations are discrescretized using fourth-order central finite 
difference scheme. At the solid boundary second-order forward differencing and at the 
outer boundary and second-order backward differencing are used. The second-order 
accurate forward difference formula at j = 1 


d<t> __ -3 d>j + 4<ftj + i - d> J+2 
dr/ 2Ar] 


d 2 <f> _ </>j — 2<p J+ i + (f>j + 2 
dr] 2 At / 2 


The second-order accurate backward difference formula at farfield ( j = N ) are 


dcp 

3 <t>j — 40j_i + 4>j - 2 

(2.53) 

dr] 

2A r] 

d 2 <p 

<f>j ~ 2<Aj-l + <j>j-2 

(2.54) 

dr] 2 

<1 


At the grid point next to the solid boundary, j = 2, equations are descretized by the 
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third-order four-point finite difference scheme 

d± _ -h 6<ftj - 3^ i+ i — 2<f> j +2 

dt] 6 At] 

d 2 d> _ <Pj — '2<f>j+i + dj+2 

dr] 2 A i] 2 

Similarly at j = N — 1, third-order four-point finite difference scheme is 


(2.55) 


(2.56) 


dd> 2 4 >j ~2 ~f Z<pj-\ — 6 (fij — 

drj 


(2.57 


d 2 d d >j — ‘2d>j + i + <pj+2 


(2.58) 


dr] 2 At / 2 

At the interior points (j = 3 — ► N — 2) the equations are descretized using fourth- 
order central difference formula 


d(f> _ — <ft/+2 + + <ftj - 2 

dr] 12At/ 


(2.59) 


d 2 ^ d*j+2 + 16(/i>j+i - 30^j + 16^-1 - <pj—2 . 6Q . 

dr/ 2 ~ 12 At/ 2 ' 1 ’ 

The descretized system of equations yields an algebraic system equations of 

the form 

AL2 <P»j_2 + + AD$j + AUl$j + j + AU2 0j +2 = 0 ; 

j = 2, N — 1 ■ (2.61) 

Here AL2, AL1, AD, AU1 and AU2 refer to the lower subdiagonal, sub- 
diagonal, diagonal, superdiagonal and upper superdiagonal matrices of size (5x5) 
respectively. <f>j represents the vector q at grid point j. The descretized system of 
equations and the homogeneous boundary conditions at the wall and at the outer 
boundary yield a homogeneous block penta-diagonal system of the equations as shown 
in figure 2.3 . 
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2.2.2 2D Eigenvalue Problem 

In this section, the numerical procedure that is used to solve the 2D eigenvalue 
problem is described. The linear stability equations for this case is given in equation 
Eq. 2.37. As it was discussed previously, the eigenfunctions are now functions of 
azimuthal [x 2 ) and normal ( 2 * 3 ) coordinates. The derivatives in the normal direction 
can be clescretized using the fourth-order central difference scheme that is used in 
the solution of ID problem. The problem is how to descretize the derivatives in the 
azimuthal direction. First, they were descretized using the central finite difference 
scheme. The solution that was obtained were very oscillatory and it was necessary 
to distribute too many grid points in the x 2 direction. This turned out to be very 
expensive and required enormous memory. Therefore, this method was abandoned 
and it was decided to use the Fourier series method to resolve the variables in the x 2 
direction. 

In this method, the dependent variables q and the coefficients of the partial 
differential equations A 2 , A 3 , A 2 3 , B 2 , B 3 and C are represented by Fourier series in 
the form 

n= oo 

q(*2,* 3 ) = E qnp 3 )e inl2 , (2.62) 

n=— oo 

771 — 00 

A(x 2 ,x 3 ) = E A m p 3 ) e ‘ mX2 1 - 7 T<.T 2 < 7 T- (2.63) 

m=— oo 

In the numerical procedure, it is necessary to replace the infinite Fourier series by 
finite sums in the form 

A/max 

q(x 2 ,x 3 )= E qn(z 3 )e‘ ni2 , (2-64) 

n= Nmax 

771 — max 

A(x 2 ,x 3 ) = E A m (x 3 )e ' mr2 , - 7 T < x 2 < 7 T , (2.65) 

771— — Xf max 

where N max and M max are the maximum number of modes kept in the Fourier expan- 
sions for the disturbances and their coefficients respectively. The Fourier components 
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A m are obtained using the Discrete Fourier Transform formula 


Am — 


1 


Mm cur “ 1 


2 A L 


E A m (x 2 ,x 3 )e 


imx2 


(2.66) 


Tn — — Mm ax 

The derivatives of the eigenfunction q(.r 2 ,x 3 ) in the x 2 direction now become 

dq(x 2 ,.r 3 ) 

dx 2 


= E ^qn(x 3 )e ,ni2 , 

Tl= — Nmax 

(2.67) 

— A' m rt x 

= E (-n 2 )qn(x 3 )e , “ I • 

71— iV max 

(2.68) 


dxS 

Substituting these expressions into the governing equation and collecting terms, one 
obtains the following ordinary differential equation for each Fourier mode, n 0 . 


E A-3( no _n) + (in A.23(n 0 -n) + B3( no _ n )J 


71 — N 2 

E 

n=N l 


d q n 

dx:\ 


+ ( — n 2 A2(n 0 -n) + i n B2(n 0 -n) + C( no _ n )| q n — 0 , ( 2 . 69 ) 

n Q — Nmax, A max- 


Here 


N \ — ■ 7TlZTl{ Nmaxi TXq T -A/j^ax} , 

N 2 = max{—N max ,n 0 — M max } ■ (2.70) 

Since the mean velocity is symmetric about the windward plane, x 2 = 0, 
the stability equations permit symmetric and antisymmetric type disturbances. For 
symmetric disturbances, the Fourier modes are related by 


{Al(_„), U 3 (_ n ), $-n) — {^l(n)i ^3(n)5 Pni @n} i 

U 2 (_ n ) - — W2(n)i n = 0, Nmax 


(2.71) 

(2.72) 
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and for antisymmetric disturbances 

{ttl(— n)i — P—n i $ —n } = f^3(n) t Pm } i (2.^3) 

^2 ( — n) = ^2(n)i ^ = 0, N max ■ (2.(4) 

Hence, it is sufficient to solve the equation, Eq. 2.69, for n 0 = 0. N max . Now, 
the equations Eq. 2.69 which is similar to equation Eq. 2.44 in ID problem and can be 
solved using the fourth-order central difference scheme. When the equation Eq. 2.69 
is descretized using the fourth-order central formula, again an algebraic system ot 
equations in the form 


AL2 $j_ 2 + AL1 $j_! + AD<f>j + AU1$ J+I 4- AU2<P J+2 = 0 ; 

j = 2, N - 1 • (2.75) 

is obtained. The size of the matrices AL2, ALl, AD, AU1 and AU2 now becomes 
{5 x{N max + l),5x(N max + 1)}. The descretized system of equations and the homo- 
geneous boundary conditions at the wall and at the outer boundary again yield a 
homogeneous block penta-diagonal system of equations as shown in figure 2.3. 

Eventhough the block penta-diagonal system can be solved efficiently, it was 
found that if the system is rewritten as a banded system it can be solved two times 
faster using the LAPACK subroutines ZGBTRF and ZGBTRS. The transformation of 
penta-diagonal matrix system into a banded matrix is easily implemented as quoted 
in LAPACK user’s guide. 

If A represents the penta-diagonal matrix system i.e, 


{AD} 

[at/ii 

[AU 2] 



[Ail] 

[AD] 

[AC1] 

[AU 2] 


[AL2\ 

[All] 

[AD] 

[AU 1] 

[AU 2] 




[AL 2] 

[All] [AD] . 
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minimum requirement of 50 Fourier modes in the azimuthal direction and 30 points 
in the wall-normal direction, the size of the matrix A becomes (5x50 )x30 = 7500 
and it takes about 32 hours of CPU-time on a Sun-Ultra 2 workstation. 

However, in real situations it is not required to find all the eigenvalues. 
It is sufficient to compute a few most unstable eigenmodes. This can easily be 
obtained using Implicitly Restarted Arnoldi Method using ARPACK software package 
developed by Lehoucq, Sorensen, Yang [lj. This method finds specified number of 
eigenvalues in a region close to a specified point in the complex u> plane very efficiently. 
When applied to an eigenvalue problem with the matrix A of a size of 10000, it takes 
only 20 minutes on the same workstation mentioned above to compute 10 eigenvalues 
located close to the specified point in the complex u; plane. 

The rest of the section will focus on the details involved in the formulation of 
generalized eigenvalue problem to be solved using the Arnoldi Method. It is observed 
from linear stability equation Eq. 2.34, that the frequency uj appears as simple ' rms 
of first power. 

- Continuity equation : — iujp + ■ • •, 

- xj momentum equation : —iupui + • • •, 

- X 2 momentum equation : —iujpU 2 + • • 

- X 3 momentum equation : — iuipu 3 + • • 

- Energy equation : —iu>p6 + • • • 

Therefore, in an attempt to capture the eigenvalues accurately, the global method 
is formulated as a temporal one. The temporal global eigenvalue formulation is as 


follows. 
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In order to solve an eigenvalue problem like equation Eq. 2.78 using the 
Arnoldi algorithm using the software package ARPACK, it is required that the ma- 
trix B be a symmetric positive definite or a symmetric semi-definite. The appearance 
of temporal frequency wasa first power terms only, in the stability equations, suggest 
an appropriate selection of disturbance vector of q = {ux,u 2 , U 3 , p, 0 } T with corre- 
sponding order of stability equations - a momentum, x 2 momentum, ,r 3 momentum 
equations, continuity equation and energy equation respectively. This assures the di- 
agonal, and hence symmetric formation ol matrix B. This is the reason for replacing 
pressure fluctuation tt by density and temperature in the stability equations. 

A further reduction of the generalized eigenvalue problem Eq. 2.78 to a 
more storage efficient standard eigenvalue formulation of the form 


A3> = u;3> , 


(2.79) 


can be obtained by manipulating the homogeneous boundary conditions. The pro- 
cedure will be explained for the case of ID eigenvalue formulation, and the similar 
manipulation can be applied to 2D eigenvalue as well. The only zero rows in matrix 
B are that corresponds to the boundary conditions Wi,u 2 ,U3 and 6 as shown below. 
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(2.80) 



And the middle elements of B take the form 


B = 


ip 

ip 

ip 


i 



(2.81) 


Without any loss of generality in the boundary conditions, the zero terms 
in B can be replaced by unity. Further, dividing the continuity equation by i and 
momentum and energy equations by ip will make B an identity matrix. Hence the 
resulting eigenvalue problem is reduced to equation Eq. 2.79. 


2.2.4 Local Method 

The local methods are used to confirm and refine the eigenvalues obtained from 
global solver. But, global methods are computationally much more expensive than 
the local methods since they compute the whole or part of the eigenvalue spectrum 
of the descretized system. However, the local eigenvalue solvers require a guess for 
the eigenvalue and hence make the use of a global solver inevitable. 

It was shown earlier that, the system of descretized stability equations can 
be formed as a lock penta-diagonal system or as banded matrix format. It can be 
seen that these equations are homogeneous and in order to avoid trivial solutions, 
Malik [34] suggests replacement of boundary condition Ui( r3= o) = 0 by the normaliz- 
ing condition that the pressure fluctuation (or equivalently density fluctuation in the 
present case ) 7 r( X 3 =0 ) = 1. Thus, the equations are transformed into an inhomoge- 
neous system. But here, a different approach will be followed. After experimenting 
with a number of different normalizations, it was found that a faster rate of con- 
vergence resulted by replacing continuity equation in the middle of the boundary 
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layer wall-normal location by the normalizing condition iq = 1 . ( In the case of 2D 
eigenvalue problem, a Fourier mode <q n = 1 is used for normalization). 

Consider the missing continuity equation denoted by »/>( A) = 0. Here A refers 
to the exact eigenvalue i.e, the streamwise wavenumber a if it is spatial formulation 
or temporal frequency u if the problem is considered as temporal. 

Suppose A 0 represents the guess for the eigenvalue and AA the error from 
the exact value so that ip( A + AA) = 0. Now, using Taylor series expansion for A A 


MK) + 


AA = - 


du' 


dAX 

D(A 0 ) 


AA = 0 , 


(2.82) 


(2.83) 


(d^/dAX) 

The iterative procedure for local ID temporal stability will be described 


below. 


* For a specified wavenumber a, formulate the penta-diagonal system for shape- 
function q, normalize the with tq = 1 by replacing continuity equation at 
wall-normal grid point j m id where the phase speed C = u} r /a T is about 0.7 . 

* Iterate on the guess value for temporal frequency u 0 until the missing continu- 
ity equation is satisfied. The correction Au is determined from the equation 
Eq. 2.83. 

The procedure for spatial formulation is similar; where temporal frequency 
ui is fixed and the Newton- Raphson iteration is done on guess value for a until con- 
vergence is reached. For 2D eigenvalue problem, a Fourier component with maximum 
amplitude is used for the normalization. 



2.3 Summary 

To conclude this chapter, the formulation are summarized as a schematic diagram. 



Fig. 2.4 Stability problem formulation 








CHAPTER 3 

MEANFLOW COMPUTATIONS 
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The stability computations basically involve two steps. The first step is to 
compute the mean velocity profiles accurately. The mean flow can be obtained by 
solving the boundary layer equations, Parabolized Navier-Stokes (PNS) equations or 
full Navier-Stokes equations depending on the problem that is being analysed. In 
this work, we investigate the linear stability of supersonic boundary layers over a 
sharp axi-symmetric cone at an angle of attack. The mean flow was computed using 
the well developed TLNS3D code which was developed at NASA Langley Research 
Center by Vatsa and Wedan (1990) [3]. 

3.1 TLNS3D 


Three-dimensional time-dependent Thin-Layer Navier-Stokes equations are used in 
TLNS3D for modelling the flow. The set of equations are obtained from the com- 
plete Reynolds-averaged Navier-Stokes equation by retaining only the viscous diffu- 
sion terms normal to the solid surfaces. The TLNS3D code employs eddy-viscosity 
hypothesis to model turbulence; the Baldwin-Lomax turbulence model is used for 
turbulence closure. In this work, the mean flow is computed for laminar flow. 

The steady-sate solutions to Thin-Layer Navier-Stokes equations are ob- 
tained using a semidiscrete cell-centered finite- volume algorithm, based on a Runge- 
Kutta time-stepping scheme. In order to suppress odd-even decoupling and oscil- 
lations in the vicinity of shock waves and stagnation points, a linear fourth or- 
der difference-based and nolinear second order difference-based dissipation is added. 
TLNS3D incorporates both the scalar and matrix forms of the artificial dissipation 
models; but in the mean-flow calculation over the sharp cone matrix form is used. 



17 



Fig. 3.1 The coordinate system for the cone located in a supersonic flow at an angle 
of attack (Moo = 3.5). 

For the sharp cone the Cartesian coordinate system (x, y, x) is located at the 
vertex and a body-oriented coordinate system (x!,x 2 ,x 3 ) fixed in time is considered 
with xj along the generator, x 2 in the azimuthal direction and x 3 normal to the 
surface as shown in figure 3.1. 

The coordinate trasformation is given by 

Xj = xcos(Q) + Rsin(Q) , 
tan(x 2 ) = y/z , 

x 3 = Rcos(Q) — xsin(Q) , 


(3.1) 

(3.2) 

(3.3) 


where 

r ?2 _ „2 , _2 

n = y + z ■ 

The govering equations can be written in the conservation form for the body-fitted 
coordinate system as 
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0 /T _ lrM dF dG 3H 9G v 

— {J U) + -r f x V -x — -X , 

at ox i dr 2 ox 3 dr 3 


(3.4) 


where U is the conserved varibale vector and F, G , and H represents the convective 
flux vectors. G v represents the viscous flux vectors normal to the body surface. Here 
only the viscous diffusion terms in the ,C3-direction are retained, due to the fact that 
in high-Reynolds-number flows, dominent contribution to the viscous effects are from 
viscous diffusion normal to the body-surface. 


u = { 


p 

pu 

pv 

pw 

pE 


F = J 


-i 


G = J~U 


pu 1 

PU\U + X\ x p 
pU\V + XlyP 

pillW + x Xz p 
pu\H 

pu 2 

pu 2 u 4- X 2x p 

PU 2 V -(- X 2y P 
pu 2 w -I- X 2z p 

pu 2 H 

PU3 

pu 3 u + X 3x p 
pn 3 v -f X 3y p 
pu 3 w + x 3z p 

pu 3 H 


(3-5) 


(3.6) 


(3.7) 


H = J -1 < 


(3.8) 
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r - 

i? eoo J 


0 

©1 + X.3 x ©2 

©I brj + X 3y (f) 2 ► f 

©1 + X 3z (p 2 

<P\d + €’©2 


(3.9) 


where 


©l — r 3i- 2 + .P.3J/ 2 + ;i'3 s 2 , 


©2 — x(-C.3a-«r 3 + X 3y V r ^ + X 3z W Xs ) , 


g 2 = u 2 4- w 2 + to 2 , 


(3.10) 


V 2 ),3 + (^- i )^ r -’ 


Here (u,t>,tt>) are the velocity components in the (x,y,z) directions and 
(ui, u 2 , U 3 ) are the contravarient velocity components (velocities along the body-fitted 
coordinate directions) which are defined as 


III = Xi x U + Xi y v + Xi z W , 
u 2 = x 2x u + x 2y u + x 2z w , 

u 3 = x 3x u + x 3y v + x 3z w , ( 3 . 11 ) 


p is the density, p is the pressure, and E is the total energy. J is the Jacobian of 
the transformation and X{j are the direction cosines. Additionally, p is the viscosity, 
A /00 is the free-stream Mach number, Re^ = -f- is the free-stream unit 

Reynolds number, and Pr = is the Prandtl number. 

Distribution of temperature can easily be obtained from the expression 


T = 


E — ~ ( u 2 + v 2 + tv 2 ) 
£ 


7(7 - 1)^00 • 


( 3 . 12 ) 



Since the steady state flow is symmetric about the plane through the wind- 
ward and leeward rays, the mean flow computations were done in the half plane ( 
,r 2 = 0° to X 2 = 180° ). By symmetry. 


{u,,u 3 ,p,r}U = {(ii, u 3 ,p,r}|_ l2 , (3.13) 

^^2 |x-2 = ^2 ' (3.14) 

3.2 Boundary- layer Profiles 

The physical parameters for which the computations were performed are given in 
table 3.1 . 


Table 3.1 Details of the cone problem 


Parameter 

Value 

Half cone angle 
Angle of attack 
Free-stream temperature 
Free-stream Mach number 
Unit Reynolds number 
Cone length 

5° 

2° 

94 °K 
3.5 

8.7xl0 6 /m 

1.574m 


The computations were done at the adiabatic wall conditions. In the com- 
putation of the meanflow using TLNS3D a mesh size of (97x257x49), i.e, 97 points 
in the streamwise direction, 257 points in the normal direction and 49 points in the 
azimuthal direction ( windward ray to leeward ray) were used. 

The mean-flow profiles at four different streamwise locations (:ri = 0.033m, 
0.197m, 0.3505m and 0.5105m), where the linear stability computations were per- 
formed, are presented. The thickness of the boundary layers at the four streamwise 
locations for different azimuthal stations are plotted in figures 3.3. Here, the bound- 
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Fig. 3.2 Schematic diagram of cross-flow and streamwise velocity profile (Reed Sz 
Saric). 

ary layer thickness is defined at a wall-normal location where the values of streamwise 
velocities at consecutive normal points differ by less than 0.1%. It is seen that the 
boundary layer thickness inceases sharply about three times when the flow goes from 
leaward to windward plane. In the left ordinate the boundary layer thickness is plot- 
ted in mm and in the right ordinate the nondimensional thickness £3 is shown, which 
is defined as 

(3.15) 

Figure 3.4 depict the variation of the boundary layer edge velocity U e in the az- 
imuthal direction at different streamwise locations. As it is seen the variation of the 
velocity is small in the azimuthal direction. 

First, the meanflow profiles at the streamwise station x x = 0.033m are 
given. The velocity profiles tangent to the inviscid stream line at different azimuthal 
locations are shown in figure 3.5. There are significantly different characteristic 
profiles on the windward and leeward planes. On the the windward side, boundary- 


£3 = 


*3 


*1* 
V Ue* 
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layer is comparatively thin and gradients are quite larger near the wall. On the 
leeward plane much thicker layer with smaller surface shear is apparant. 

Figure 3.6 show the cross-flow profiles along azimuthal x 2 locations. The 
cross-flow components were non-dimensionalized with the boundary- layer edge ve- 
locity, and positive cross-flow is taken towards positive x 2 direction and awav from 
zi direction as shown in figure 3.2. The cross-flow velocities exhibit local maximum 
values at x 2 = 90° and reaches values between 4% and 6% of the freestream velocity. 
Closer to the leeward side (around x 2 = 160° ) negative cross-flow is evident i.e, 
cross-flow component changes in sign before it decays to zero towards the boundary- 
layer edge. As could be expected, cross-flow phenomena dissapear in windward and 
leeward rays due to the symmetry in the mean-flow. 

Figure 3.7 shows the azimuthal velocity distributions at different azimuthal 
locations. Figures 3.8 and 3.9 show the density and temperature distributions at 
this x i location. The similar results of the mean-flow profiles at the other three 
streamwise locations z x = 0.3505m and 0.5105 m are presented in figures 3.11 to 
3.19. 

The velocities are nondimensionalized by the characteristic velocity U* e 
which is taken to be the local velocity at the edge of the boundary layer in the 
streamwise direction in ID formulation and the boundary layer edge velocity at a 
reference station x 2 — 90° in the case of 2D formulation. The normal coordinate is 
normalized by the length scale 

(3.16) 

where v* e is the kinematic viscosity at the edge of the boundary layer. It is seen 
that except close to the leeward plane region the boundary layer profiles are almost 
linear in most part of the boundary layer and the boundary layer thickness increases 
gradually. Very close to the leeward plane , the profiles exhibit very strong inflectional 
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character, especially 165° < x 2 < 180°. Figure 3.6 shows the crossflow velocity 
component in the positive x 2 direction at different azimuthal locations. It is observed 
as expected that the crossflow velocity is zero at the windward and leeward planes 
and reaches maximum in the middle of the refion 90° < x 2 < 120° The maximum 
crossflow velocity is about 4% of the boundary layer edge velocity. From this it is 
also expected that the boundary layer will be very unstable in the middle region 
compared to the region near the windward or leeward planes. 

Figures 3.21 shows the contours of the crossflow Reynolds number defined 


bv 


Re r = 


Qmax&\0% 


(3.17) 


in the x x x 2 plane. Here Q max is the maximum cross-flow velocity (located at 6 max 
), and 8 is the thickness defined by the point above 8 max at which the cross-flow 
velocity is 10% of Q ma x- As expected the cross flow Reynolds number is maximum 
in the middle region and increases with the x 2 direction. The maximum cross-flow 
Reynolds number is about 250 at £j = 0.5105m. 



150 


a) x\ = 0.033m 


b) x\ = 0.1978m 




a) xi = 0.3505m b) = 0.5105m 


Fig. 3.3 The variation of boundary-layer thickness in the azimuthal direction at 
different streamwise locations. 







a) x\ = 0.033m 


b) x\ = 0.1978m 




a) xi = 0.3505m b) x x = 0.5105m 


Fig. 3.4 The variation of boundary-layer edge velocity in the azimuthal direction at 
different streamwise locations. 
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Fig. 3.10 The distribution of Fourier components of a) streamwise velocity, b ) az- 
imuthal velocity and c) temperature with Fourier modes m. 
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Fig. 3.17 Cross-flow velocity profiles. (Re = 2211, £i — 0.5105m). 
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CHAPTER 4 

RESULTS AND DISCUSSION 

In this chapter, the linear stability results obtained using the ID eigenvalue 
method and the 2D eigenvalue method are presented. In order to study the evolution 
of the disturbances downstream, the ID and 2D stability computations are performed 
at different locations along the streamwise direction : X\ = 0.033m, 0.197m, 0.3505m 
and 0.5105m. The results are presented in two sections. In the first section, the 
computational results from ID eigenvalue approach are given and the results from 
the 2D eigenvalue approach are presented in section 2. 


4.1 ID Eigenvalue problem 

Before proceeding to present the linear ID eigenvalue results, the computational grid 
used in the computations are described. As mentioned in Chapter 2, the height of the 
computational domain was taken to be about four times the boundary layer thickness. 
This was necessary because, eventhough most of the perturbations decay to zero 
within the boundary layer, density and normal velocity perturbations persist until 
about three to four times the boundary layer thickness. The computational domain 
consisted of 85 points in the wall-normal direction, with first 43 points clustered 
within the boundary layer according to the algebraic grid-stretching given by equation 
Eq. 2.47. Because of enormous requirement of the storage, the grids points in the 
2D eigenvalue computation are limited to 49 points. To create a general platform for 
the comparison of linear ID and 2D eigenvalue results, ID eigenvalue problem was 
also solved on a grid size of 49 for few cases. The stability results computed on a 
grid size of 49 are found to be accurate within five decimal points to that of obtained 
with 85 grid points. 



It was explained in Chapter 2 that in the eigenvalue computations, since 
the local method requires a guess value, a global method has to be first used to 
compute the whole or part of the eigenvalue spectrum. In the present work, the ID 
global eigenvalue computations are performed using the QZ algorithm of the subrou- 
tine ZGEGV in LAPACK software library. For a typical ID problem with 85 points 
in the normal direction, the leading dimension of the matrix A of the generalized 
eigenvalue problem becomes (5x85) = 425 and it takes only 30 seconds on a Sparc- 
Ultra-2 workstation (333MHz) to compute all of the eigenvalues. These spectrum 
of eigenvalues showed only one, or two physical eigenvalues which are unstable. It 
was found that, in ID eigenvalue computations, the initial guess value need not be 
close to the unstable eigenvalue sought. Therefore, ID eigenvalue computations are 
performed, in most cases, using the iterative local solver by the continuation method, 
i.e, using the eigenvalue at the previous location as the initial guess for the current 
location. 

To orient the reader with the co-ordinate system used, a schematic diagram 
showing flow and wave propagation directions is depicted in figure 4.1. Here m which 



Fig. 4.1 Schematic diagram of the flow and the wave directions 
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is an integer, is the number of waves in the azimuthal direction. In order to determine 
the unstable region of the boundary layer, first several eigenvalue computations were 
performed along the azimuthal direction from windward ray to the leeward ray at 
number of fixed streamwise locations aq, for different streamwise wavenumber a and 
azimuthal mode number m. A sequence of such computations starting from a aq 
location closer to the tip of the cone exhibited that the region up to aq = 0.03m 
is linearly stable. Hence, the linear stability results at four aq locations starting 
from aq = 0.033m are given here. The length, velocity and time scales used to 
nondimensionalize the variables are as follows. 

Table 4.1 Parameters used in nondimensionalization of different scales. 


Parameter 

Streamwise distance, X\ /[m\ 

0.033 

0.197 

0.3505 

0.5105 

L /[mm] 

0.0627 

0.146 

0.192 

0.2.31 

U e /[m/s] 

663.3 

659.7 

658.4 

658.2 

tx 10 8 /[s] 

9.46 

22.13 

29.18 

35.11 


The linear ID stability results at streamwise locations aq = 0.033m, aq = 
0.197m, aq = 0.3505m and aq = 0.5105m are shown in the figures 4.2 to 4.6, 4.8 
to 4.11, 4.12 to 4.19 and 4.20 to 4.22 respectively. Because the stability results 
at these stations are similar, only the results at aq = 0.033m and aq = 0.3505m are 
discussed in details. 
aq = 0.033m 

The figures 4.2a, 4.2 6, 4.3 a and 4.36 show the variation of temporal ampli- 
fication rates u t and temporal frequency u T with streamwise wavenumber a for dif- 
ferent azimuthal mode number m computed at the streamwise location aq = 0.033m 
at different azimuthal positions aq = 45°, aq = 90°, aq = 120° and aq = 150° 
respectively. The results show that the temporal frequency uj t varies linearly with 




wavenumber a for all the azimuthal mode numbers shown and that the phase speeds, 
defined as C = are approximately constant and equal to 0.7. Further, the station- 
ary waves, i.e, ay = 0, are stable at this ,Ci location. At .r 2 = 45°, the disturbances 
are linearly unstable for —15 < m < -9. The maximum amplification rate is about 
0.0006 and occurs at m = -12 and a = 0.1 . The similar variations of the temporal 
amplification rates are observed at x 2 = 90°, x 2 = 120° and .r 2 = 150°. However, 
as one moves from the windward side .r 2 — 45°, toward the leeward side x 2 = 150°, 
the unstable streamwise wavenumber a shift from a range of (0.07 ~ .13) to (0.02 ~ 
0.07). The locally most amplified frequency is found to be a’ = (0.038,0.00101) at the 
azimuthal location of x 2 = 120° for a — 0.0c and m — 9 and lo — (0.03,0.00103) 

at x 2 = 150° for a = 0.05 and m = -6. The negative sign of the mode number m 
indicates that the propagation of the disturbance wave is in the negative side of x 2 . 

The local wave propagation direction is given by 

ip ■ tan 1 — (4-1) 

Oc r 

where R is the radius of the cone. The wave angles for the most amplified wave at 
different x 2 locations are given in the Table 4.2. In this table, e is the inclination of 
the inviscid stream line to the aq axis (figure 4.1). One can see that in the negative 
x 2 direction most amplified waves are inclined at about 70° from the inviscid stream 
line. Figures 4.4 and 4.5 show the variation of the amplification rate u>i and the 

Table 4.2 The waveangles for the most amplified wave at different x 2 locations. 


*2 /[°] 

e/[°] 

-0/H 

<f> = —ip -|- e 

45 

1.68 

69.09 

70.77 

90 

2.73 

69.09 

71.82 

120 

2.30 

71.67 

73.39 

150 

1.22 

71.02 

72.24 


frequency ay along the azimuthal direction for a constant wavenumber a — 0.05 
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for different azimuthal mode number m. Figures 4.6 and 4.7 depict the results 
for a = 0.07. It is observed that the frequencies are lower in the middle than they 
are near the windward and leeward regions. This is due to the fact that cross-flow 
instability is dominated by low frequency disturbances compared to that without 
the cross-flow. As discussed in Chapter 2, the cross-flow is maximum in the middle 
region and hence the most unstable disturbances are low frequencies compared to 
that near the windward and leeward plane. 
gi = 0.3505m 

The variation of the temporal amplification rate u.\ and frequency u r with 
wavenumber a r at x 2 = 0°, 20°, 41°, 97°, 120° and 160° at the streamwise location 
X\ = 0.3505m are plotted in figures 4.12a, 4. 126, 4.13a, 4.136, 4.14a and 4.146. 
Figures 4.15 and 4.16 show the variation of the amplification rate and the frequency 
along the azimuthal direction for a constant wavenumber a = 0.07 for different 
azimuthal mode number m. It is seen from the figure 4.12 that the amplification rates 
are the highest closer to a — 0.07. First observation is that the amplification rates are 
high for disturbances with m = —30 to —40. The maximum amplification rate is = 
0.0046 and this occurs around X 2 = 90° for m = —30 and a = 0.07. The amplification 
rates between 40° < x 2 < 130° vary in the range from 0.004 to 0.0046 and it decreases 
gradually to 0.0026 and 0.003 at x 2 = 0° and 160° respectively. The most unstable 
frequencies are lower in the middle region compared to that near the windward and 
leeward planes. For m = —30 and a = 0.07 the frequency of unstable disturbance 
is 0.039 at x 2 = 80° and they are 0.046 and 0.051 at x 2 = 0° and 160° respectively. 
These translate to 17.67 kHz at x 2 = 80° and 20.83 kHz and 23.1 kHz at x 2 = 0° and 
160° respectively. It is also seen that the frequencies of the unstable disturbances 
decrease with increasing negative m. For ra = —60, the unstable frequency is 0.034 
at x 2 = 80° and it is 0.041 for m = -20. The reason for this is that with increasing 
m, the wavevector aligns closer to the cross-flow direction and the frequencies of the 
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unstable disturbances decrease. The results show that the amplification rate of the 
disturbances with positive m. are small except closer to the windward plane. It is 
also observed that at the windward plane x 2 = 0° the eigenvalues for the positive 
and negative m have the same values. This is due to the fact that at x 2 = 0°, 
the cross-flow velocity is zero and the meanflow is two-dimensional and it does not 
differentiate between positive and negative m values. When one moves away from the 
windward plane the cross-flow velocity increases, hence the mean flow becomes three- 
dimensional and the eigenvalues depend on the positive or the negative direction at 
which the wavevector is aligned. Another observation is that there is an apparent 
symmetry about the eigenvalues about x 2 = 90°. This can be explained from the 
observation that the maximum cross-flow velocity increases from zero to a peak value 
in the middle and decreases to zero again at the leeward plane. Hence the velocity 
profiles are approximately symmetrical about the middle plane and it is expected 
that the eigenvalues will also be symmetric. This observation becomes important 
when the results from 2D eigenvalue approach is interpreted in the next section. 

As a prelude to later comparisons, all the eigenvalues obtained for different 
azimuthal mode numbers m at different azimuthal locations for the wavenumber 
or = 0.074 are plotted in the complex u-plane in figure 4.17. This is a representation 
of figures 4.15 and 4.16 in the complex plane. Some of the azimuthal locations x 2 at 
which these eigenvalues are computed are also marked. One thing to conclude from 
this figure is that the eigenvalues are clustered in a confined region in the complex em- 
plane. Figure 4.18 and 4.19 depict the amplitudes of eigenfunction distribution for 
the streamwise velocity |ui| and normal velocity |ua| at different azimuthal locations 
for m = —30 and a = 0.07. It is seen that the eigenfunctions for the streamwise 
velocity decrease to zero at the edge of the boundary layer, however, the eigenfunction 
for the normal velocity decrease to zero slowly. Another observation is that the 
locations of the maximum amplitude move towards the edge of the boundary layer 
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when one moves toward the leeward plane from the windward plane. Similarly fig- 
ures 4.20 to 4.22 show the ID stability results obtained at the streamwise location 
xi = 0.5105m. 

Summarizing the linear stability results of ID eigenvalue method the fol- 
lowing conclusions can be made. 

• The boundary layer region up to ;t x = 0.03m is linearly stable and is unstable 

beyond that location. However, the neutral stability region is not a straight 
line across the cone at art = 0.03m but is curved with the front of the neutral 
stability region falling near 90°. This is due to the effect of varying cross-flow 
components from windward to leeward locations. 

• The most amplified temporal amplification rate occurs around x 2 = 90° for a 

streamwise wavenumber a = 0.07 and azimuthal mode number m = —30, —40. 

• The effect of cross-flow component is dominant in the middle region in azimuthal 

direction and this is manifested in the increase in the temporal amplification 
rate around x 2 = 90° for negative m — — 10 to —60 and the temporal amplifi- 
cation rate for positive m values. 

• The unstable temporal disturbance waves with most amplification rate u>, travel 

in the negative side of x 2 direction. The unstable waves propagating in the 
positive side of x 2 become stable after x 2 = 90°. 




Fig. 4.2 :The variation of temporal amplification rate u) t and temporal frequency 
with streamwise wavenumber a. (xi = 0.033m). 
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Fig. 4.3 :The variation of temporal amplification rate and temporal frequency u 
with streamwise wavenumber a. (aq — 0.033m). 
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Fig. 4.4 : Distribution of temporal amplification rate (u>i) with x<i (streamwise 
location x\ = 0.033m, a = 0.05). 
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Fig. 4.6 :Distribution of temporal amplification rate u;,- with X 2 ■ ( 2:1 = 0.033m, 
a = 0.07). 
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Fig. 4.10 Distribution of temporal amplification rate (u>i) with x-i- (stream- 
wise location x\ = 0.197m, a = 0.07) 
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b) x 2 = 20° 


Fig. 4.12 The variation of temporal amplification rate and temporal frequency u r 
with streamwise wavenumber a. (xi = 0.3505m). 
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Wavenumber a r 


a) x 2 = 41° 



Fig. 4.13 The variation of temporal amplification rate a;,- and temporal frequency w, 
with streamwise wavenumber a. (xi = 0.3505m). 
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b) z 2 = 160° 

Fig. 4.14 The variation of temporal amplification rate u and temporal frequency t J, 
with streamwise wavenumber a. (xi = 0.3505m). 
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Fig. 4.15 Distribution of temporal amplification rate (u;,-) with x 2 . (streamwise 
location x x = 0.3505m, a = 0.07). 
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Fig. 4.20 The variation of temporal amplification rate 1 and temporal frequency to, 
with streamwise wavenumber a. (x x = 0.5105m). 




Fig. 4.21 :Distribution of temporal amplification rate u?,- with x 2 . 
(xi = 0.5105m, a = 0.07). 



Fig. 4.22 :Distribution of temporal frequency (u> r ) with x 2 (streamwise location 
xi = 0.5105m, a = 0.07). 
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4.2 2D Eigenvalue Problem 


In this section, the results from the 2D eigenvalue computations are presented. The 
computations are performed at different aq locations xj = 0.033m, Xj = 0.197m , 
Xj = 0.3505m and x\ = 0.5105m. 

As mentioned earlier, the main difficulty in the 2D eigenvalue computations 
is the requirement of large memory capacity. Therefore, in a typical 2D eigenvalue 
computations the maximum number of grid points in the wall-normal direction is 
limited to 49 points and the Fourier modes to a maximum of 89. However, in the 
regions close to the tip of the cone, it was found that most of the unstable disturbances 
could be captured accurately with 59 Fourier modes. This permitted the eigenvalue 
computations at aq = 0.033m be performed with 65 wall-normal points and these 
results are found to be consistent with that of 49 normal points and 89 Fourier modes. 

Unlike in ID stability computations where the eigenvalue spectrum for a 
given azimuthal mode number m consisted only of a few sparse unstable eigenval- 
ues, it was found from the 2D eigenvalue computations that the eigenvalue spectrum 
showed a clustered nature of the unstable eigenvalues in the complex u> plane. There- 
fore, it is necessary that prescribed initial guess for the local 2D solver be accurate, 
otherwise the solution might converge to some other eigenvalue. For the above men- 
tioned typical 2D problem with 49 points in the wall-normal direction and 89 Fourier 
modes, the leading dimension of the matrix A of the generalized eigenvalue problem 
becomes (5x[89+l]x49) = 22050. It takes more than 200 CPU-hours on a Sparc- 
Ultra-2 workstation(333MHz) and requires about 450 MW memory to compute all 
the eigenvalues. Therefore, the ARPACK software package employing the Implicitly 
Restarted Arnold! Method (refer Chapter 2 for details) is used to obtain a specified 
number of eigenvalues in a region close to a given point in the complex u plane. 
When applied to the same problem of size 22050 , it takes only about 2 hours to 
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compute 10 eigenvalues that are located close to specified region of interest, on the 
workstation quoted above. 

The 2D stability computations performed at streamwise locations Ci = 
0.033m, 0.197m, 0.3505m and 0.5105m and the corresponding results are presented in 
figures 4.43 to 4.61 , 4.62 to 4.76, 4.23 to 4.41 and 4.77 to 4.90 respectively. From 
the linear stability results of ID eigenvalue method it was noted that the maximum 
temporal amplification rates in the unstable boundary layer region occurred for the 
range of wave numbers a = 0.06 to 0.1 . Therefore the 2D eigenvalue computations 
are performed for the same range of a. The stability results at the above mentioned 
stations are found to be similar, and hence results at = 0.3505m are discussed in 
details first. 
xi = 0.3505m 

Figures 4.23 to 4.41 show the results obtained at station aq = 0.3505m. for 
seven different eigenvalues, both symmetric and anti-symmetric modes, for a 0.0 ( . 
As it was discussed earlier, there exists large number of eigenvalues and the results 
are presented only for the most amplified disturbances. In the symmetric mode, 
the disturbances are assumed to be symmetric about the windward plane ( x 2 0 ) 

and in the antisymmetric mode the disturbances are taken to be antisymmetric ( for 
detailed explanations refer Chapter 2). 

Figure 4.23a shows the distribution of streamwise velocity disturbance Uueai 
along the azimuthal direction at the wall-normal height x 3 = 0.721mm where it 
has the maximum amplitude for the wavenumber a = 0.07 and the eigenvalue is 
uo = (0.0386,0.00437). The local Reynolds number is 1823 which is defined as 

(4.2) 

It can be seen that the disturbances are symmetric about x 2 = 0° and confined 
between 80° to 140° with the maximum amplitude occurring at 120°. The corre- 
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sponding distribution of Fourier components of the streamwise velocity disturbance 
jt/jy I with Fourier modes m is plotted in figure 4.236. It is noted that the Fourier 
components |tii/| falls in a bell-shape distribution with negligible magnitude until 
m = 20 then sharply increasing to a maximum around m = 40 ~ 50 and thereafter 
decaying gradually to zero around m = 80. The contour plots of |ui| and Ui rf , a i in the 
positive x 2 x 3 plane for the same temporal eigenvalue u are depicted in figure 4.24a 
and 4.246. The distribution of the streamwise velocity disturbance profile and nor- 
mal velocity disturbance profile are plotted in figures 4.25a and 4.256 respectively. 
It can be observed that the normal velocity perturbation u 3 persist until about four 
times the boundary layer thickness, whereas the streamwise velocity perturbation ad 
decays to zero within the boundary layer, which is equal to x 3 = 20 in nondimensional 
quantity. (The variation of the boundary-layer thickness in azimuthal direction at 
different streamwise locations are presented in figure 3.3). 

The figures 4.26 to 4.27, 4.28 to 4.29 present the results of symmetric 
disturbances with w = (0.04205, 0.00272) and w = (0.0491, 0.00264) which are similar 
to the previous results corresponding to w = (0.0386,0.00437). However, a closer 
examination of these results show that the clustered disturbances shift towards the 
leeward side with disturbance amplitude peaking at x 2 = 140° and x 2 — 160° for 
u 3 = (0.04205,0.00272) and u = (0.0491,0.00264) respectively. 

The 2D stability results for u = (0.04199, 0.00349) are shown in 4.30 to 4.32 
and, in contrast to the previous results discussed, they have some different interesting 
features. Figure 4.30a shows the distribution of streamwise velocity disturbance 
U\ rea i along the azimuthal direction at the wall-normal height x 3 = 0.487mm. An 
important observation is that , considering only the half azimuthal plane X 2 — 0 to 
x 2 = 180°, the disturbances are clustered in two different regions - around x 2 = 30° 
and 140°. This complies with the observation of ID results that there is an apparent 
symmetry of the eigenvalues about x 2 = 90° caused by an approximately symmetric 
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distribution of cross-flow in the middle azimuthal direction. However, looking at the 
figure 4.306 which depicts the distribution of u^eai along the azimuthal direction at 
.r 3 = 0.85mm, it is apparent that the disturbances clustered around .r 2 = 30° decay to 
zero whereas the counterparts at x 2 = 140° grow to a maximum. The corresponding 
distribution of Fourier components of the streamwise velocity disturbance |tti/| with 
Fourier modes m is plotted in figure 4.31. The similar results of an eigenvalue uj = 
(0.0451,0.00272) that exhibits clustered disturbance eigenfunctions around .r 2 = 10° 
and 150° are shown in figures 4.33 to 4.35. 

Next, the results for anti-symmetric modes are presented. The figures 4.36 
to 4.37 show the results of u> = (0.0437,0.0025). From figure 4.36a one can note 
that the disturbances are anti-symmetric and clustered around x 2 = —150° and 150°. 
Also, the results for the case of u> = (0.0435, 0.00315) are presented in figures 4.38 to 
4.41. The results show that the disturbances are clustered around x 2 = 30° and 150°, 
considering only the positive half of the azimuthal plane. However, the disturbances 
around x 2 = 30° are dominant and larger than the disturbances closer to the leeward 
side. 

An important observation about the results of symmetric and anti-symmetric 
disturbances is that for identical eigenvalues , there are apparently not much differ- 
ence between the eigenfunction distributions of symmetric and anti-symmetric modes 
which are clustered around the middle of the azimuthal plane. However, eigenfunc- 
tions that peak close to the windward side exhibit a noticeable differences between 
symmetric and anti-symmetric modes. This observation suggest that these clustered 
disturbances in the middle of the azimuthal direction are not affected by the mean- 
flow quantities (especially the azimuthal velocity components ) away from them and 
there is not much interaction between the disturbances in the positive half and the 
negative half of the x 2 x 3 plane. 

Finally, the distribution of the 2D eigenvalue spectrum need to be discussed. 
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Figure 4.42a shows the plot of temporal amplification rate with frequency uj t 
for wavenumber a = 0.07 at = 0.3505m. The maximum amplification rate is 
LJi = 0.00437 and the corresponding frequency u; r = 0.0386. The eigenfunctions for 
this eigenvalue are clustered in the middle of the azimuthal plane and peak around 
.r 2 = 120°. The x- 2 locations where the eigenfunctions peak are also marked on the 
figure. An eigenvalue marked with two angles indicates that the disturbances are 
clustered at two isolated x 2 regions and peak at the angles quoted. In figure 4.426 , 
the ID eigenvalues along with the 2D eigenvalues are plotted on complex u ; plane. A 
closer examination of this figure show that the most amplified temporal amplification 
rates of 2D stability method occur around .r 2 = 120° whereas in ID method they fall 
around x 2 = 90°. The shift in the most amplified temporal amplification rates u>, of 
the 2D eigenvalues toward the leeward side can be possibly explained as follows. The 
instability is determined by the degree of inflection of the meanflow profiles and the 
contribution of cross-flow components towards instability is maximum at .r- 2 = 90°. 
This is clearly manifested in the maximum increase of the amplification rates of 1 D 
eigenvalue method. However, in the case of 2D eigenvalue method which incorporates 
the variation of meanflow in the azimuthal direction, the meanflow is more unstable 
towards the leeward side and therefore shifting the most amplification rates toward 
the leeward side. 
x 1 = 0.033m 

Figures 4.43 to 4.61 depict the results obtained at x x = 0.033m for different 
axial wavenumbers a = 0.07, 0.08, 0.09 and 0.1. Figure 4.43 shows the distribution 
of eigenfunction and the Fourier components for u = (0.0391, 0.00128) and axial 
wavenumber a = 0.07. The corresponding frequency is 65.8 kHz. It can be noted 
that the eigenfunction clustered around 120°. The maximum amplification of the 
eigenfunction occurs at a wall-normal height of x 3 = 0.255m, and which is about 70 
% of the boundary-layer thickness ( For detailed data on distribution of boundary 
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layer thickness refer to figure 3.3 ). The contour plot for the absolute value of the 
streamwise perturbation velocity is shown in figure 4.44. Figure 4.43 shows the 
eigenfunction distribution for the axial velocity and the normal velocity at different 
azimuthal stations. The similar plots for an eigenvalue lj = (0.0387, 0.00311) are 
given in figures 4.46 to 4.47. In figure 4.48 the plots of eigenfunction distributions 
are given for the frequency of / = HOAkHz. Figures 4.49 to 4.50 show the results 
for the axial wavenumber a = 0.1 and u;=(0.0656,0.0001 14). Unlike the previous 
two cases, this eigenfunctions is clustered around x 2 = 0° for < 0.105mm and 
around x 2 = 130° for :r 3 > 0.105mm. Figure 4.49 shows the spectral distribution 
of the axial velocity at different heights with the mode number m. The distribution 
of streamwise velocity perturbation and temperature perturbation profiles along the 
azimuthal direction are depicted in figure 4.50 and 4.51 respectively. Figures 4.52 to 
4.53 show the similar results for the case of anti-symmetric mode with frequency of 
107.6 kHz, u> = (0.0639, 0.00094) and a = 0.1. Also the figures 4.54 to 4.57 depict the 
similar distributions for a = 0.08, u >= (0.05039, 0.00164) where the eigenfunction has 
maximum amplitude around x 2 = 0° and 140°. Similarly figures 4.58 to 4.61 show 
the eigenfunction distributions and the mode shapes for anti-symmetric disturbances 
at a = 0.08. 

,ri = 0.1978m 

Similarly figures 4.62 to 4.76 show the results obtained at x x = 0.1978m 
at different axial wavenumbers. Figures 4.62 to 4.63 depict the case of a symmetric 
mode where the eigenfunctions are clustered around 120°. In figures 4.64 to 4.67 
the distribution of eigenfunctions for a symmetric mode with frequency 42.2 kHz 
, u) = (0.0589, 0.00301) and wavenumber a = 0.09 are plotted. It can be noted 
that the eigenfunctions clustered around x 2 — 0° and x 2 = 150°. However, the 
amplitude of the eigenfunctions at x 2 = 0° are dominant for x$ = 0.287mm and die 
out gradually as one moves towards the edge of the boundary layer. Figures 4.68 
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and 4.69 show the results for u> = (0.0572, 0.00348) and figures 4.70 to 4.71, 4.72 to 
4.74 and 4.75 to 4.76 show the distribution of the axial velocity fluctuations along 
the azimuthal direction and the spectral distributions for different uj = (0.06056. 
0.00328), (0.0563. 0.00131), (0.0574. 0.00115) and (0.0546, 0.00161). It is observed 
that the eigenfunctions are clustered near the windward side for these eigenvalues. 

£i = 0.5105m 

The 2D stability results at the station x x = 0.5105m are given in figures 
4.77 to 4.90. Figures 4.77 to 4.78, 4.79 to 4.80, 4.81 to 4.82, 4.83 to 4.85. 
show the eigenfunctions and the spectral distribution for four different eigenvalues u> 
= (0.0385, 0.00445), (0.0423, 0.00351), (0.0419, 0.00439), (0.0507, 0.00288) for sym- 
metric modes. Eigenfunction corresponding to the first eigenvalue peaks around 120° 
and the spectral distribution has a Gausian shape as observed in other stations. The 
eigenfunctions corresponding to the second eigenvalue are confined to two isolated 
regions, one near the windward side between 0° - 50° and the other near the lee- 
ward side between 120° - 150°. The corresponding spectral shape shows two different 
types of distributions at the edge of the boundary layer the shape has a Gausian 
distribution and near the wall a modulated shape is observed. This observation is 
made in other stations when the eigenfunctions are clustered in two regions. Similar 
observation is made for the fourth eigenvalue. Eigenfunction for the third eigenvalue 
peaks around 20° and the spectral distribution has a Gausian shape with a long tail 
at higher mode numbers. 

Figure 4.86 to 4.88, 4.89 to 4.90 show the eigenfunctions and the spectral 
distributions for two different eigenvalues uj = (0.0486, 0.0427), (0.0517, 0.00411) for 
anti-symmetric modes. The eigenfunctions are confined to two different regions, one 
close to the windward side and the other close to the leeward side and as it was 
observed earlier, the spectral distributions show two different shapes. 

The results of the 2D linear stability results can be summarized as follows. 
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• Spectrum of 2D eigenvalues exhibits clustered nature of eigenvalues in the u> plane. 

• As expected earlier, the distributions of the disturbances along the azimuthal di- 

rection X 2 are clustered in confined regions. In the cases where the disturbances 
are clustered in two different regions in *2 direction, they are nearly symmetric 
about X 2 = 90°. 

• The most amplified temporal amplification rates occur around ,r 2 = 120°. 

• The differences between eigenfunctions of symmetric and anti-symmetric modes 

are pronounced when the eigenfunction distributions are confined near the 
windward side (x 2 = 0°). However, the differences are insignificant when the 
eigenfunctions are confined in the middle region. 
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wall-normal point x 3 = 0.721mm, where it has maximum amplitude. 



b) The distribution of Fourier components of streamwise velocity disturbance 
|m/| with Fourier modes m at different heights x 3 . 


Fig. 4.23 Case 1: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 21 kHz. ( u = (0.0386, 0.00437) , a = 0.07. ( 
Xl = 0.3505m and Re = 1823). 
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Fig. 4.24 Case 1: Contour plot of streamwise velocity disturbance, (frequency / — 21 
kHz , w = (0.0386, 0.00437) , a = 0.07 , xi = 0.3505m and Re = 1823 ). 




b) |u 3 | 

Fig 4 25 Case 1 : The distribution of disturbance profiles along the azimuthal 
direction, (frequency / = 21 kHz , u> = (0.0386, 0.00437) and a = 0.07, 1, = 0.3505m 
and Re = 1823) 
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a) The distribution velocity disturbance Ui rea / in the azimuthal direction at 
wall-normal point x 3 = 0.931mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
|ux/| with Fourier modes m at different heights x 3 . 

Fig. 4.26 Case 2: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 22.9 kHz. (cj = (0.04205, 0.00272) , a = 0.07 , 
x x = 0.3505m and Re = 1823) 





Fig. 4.27 Case 2: Contour plot of streamwise velocity disturbance |«i|. (frequency 
/ = 22.9 kHz , u) — (0.04205, 0.00272) , a = 0.07 , xi = 0.3505m and Re = 1823 ) 
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a) The distribution velocity disturbance u lrea i in the azimuthal direction at 
wall- normal point x 3 = 1.091mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
|ui/| with Fourier modes m at different heights x 3 . 

Fig. 4.28 Case 3: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 26.8 kHz. (u = (0.0491, 0.0026) , a - 0.07, 
Xl = 0.3505m and Re = 1823) 






b) x 3 = 0.85mm 


Fig. 4.30 Case 4: The distribution of eigenfunction for a symmetric mode with 
frequency / = 22.87 kHz , u = (0.0419, 0.00349) , a = 0.07 , = 0.3505m and Re 
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m 

Fig. 4.31 Case 4: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency / = 22.87 
kHz , u = (0.0419, 0.00349) , a = 0.07 , = 0.3505m and Re = 1823 ). 
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Fig. 4.32 Case 4: Contour plot of streamwise velocity disturbance, (frequency 
/ = 22.87 kHz , u = (0.0419, 0.00349) , a = 0.07 , x x = 0.3505m and Re = 
1823 ) 



b) x 3 = 0.85mm 


Fig. 4.33 Case 5: The distribution of eigenfunction for a 
frequency / = 24.6 kHz , uj = (0.0451, 0.00272) ,a = 0.07, X\ 


symmetric mode with 
= 0.3505m and Re = 
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Fig. 4.34 Case 5: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency / = 24.6 
kHz , w = (0.0451, 0.00272) , a = 0.07 , = 0.3505m and Re = 1823 ). 




Fig. 4.35 Case 5: Contour plot of streamwise velocity disturbance, (frequency 
/ = 24.6 kHz , u/ = (0.0451, 0.00272) , a = 0.07 , x x = 0.3505m and Re = 
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wall-normal point x z = 0.931mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
| Ml/ | with Fourier modes m at different heights x 3 . 

Fig. 4.36 Case 6: The distribution of eigenfunction and Fourier components for an 
anti-symmetric mode with frequency / = 23.81 kHz , u = (0.0437, 0.0025) , a - 0.07, 
x x = 0.3505m and Re = 1823) 





Fig. 4.37 Case 6: Contour plot of streamwise velocity disturbance |ui|. (frequency 
/ = 23.81 kHz , u = (0.0437, 0.0025) , a = 0.07 , x x = 0.3505m and Re = 1823 ) 





Fig. 4.38 Case 7: The distribution of eigenfunction for a symmetric mode with 
frequency / = 23.71 kHz. (u> = (0.0434, 0.00316) , a = 0.07, Xl = 0.3505m and Re 
3 





m 

b) x 3 = 0.85mm 


Fig. 4.39 Case 7: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency / = 23.71 
kHz , w = (0.0434, 0.00316) , a = 0.07 , x x = 0.3505m and Re = 1823 ) 
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Fig 4 41 Case 7: The distribution of streamwise velocity perturbation profiles along 
the azimuthal direction, (frequency / = 23.71 kHz , w = (0.0434, 0.00316) , a = 0.07, 
Xi = 0.3505m and Re — 1823) 
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a) 2D eigenvalue spectrum 
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Fig. 4.42 Case 7: Unstable eigenvalue spectrums of a) 2D eigenvalue method and 
b) ID (hollow symbols) and 2D (solid symbols) eigenvalue methods. ( a = 0.07 , 
x\ = 0.3505m and Re = 1823 ) 
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b) ID and 2D eigenvalue spectrum 
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wall-normal point x 3 = 0.255mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
\u lf \ with Fourier modes m at different heights x 3 . 


Fig 4 43 Case 1: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 65.8 kHz . (or = (0.0391, 0.00128) and or = 0.07 
, X\ = 0.033m, Re = 531). 
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Fig. 4.44 Case 1: Contour plot of streamwise velocity disturbance |ui|. (frequency 
/ = 65.8 kHz , u = (0.0391, 0.00128) , a = 0.07 , x l = 0.033m and Re = 531 ) 
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a) The distribution velocity disturbance Uireai in the azimuthal direction at 
wall-normal point X 3 = 0.215mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
|m/| with Fourier modes m at different heights X 3 . 


Fig. 4.46 Case 2: The distribution of eigenfunction and Fourier components for a+ 
symmetric mode with frequency / = 65.1 kHz . ( u> = (0.0388, 0.00311) and a = 0.07 
, X! = 0.033m, Re = 531) 




Fig. 4.47 Case 2: Contour plot of streamwise velocity disturbance |ui|. (frequency 
/ = 65.1 kHz , w = (0.0388, 0.00311) , a = 0.07 , = 0.033m and Re = 531 ) 


100 


150 180 


b) x$ = 0.276mm 


Fig. 4.48 Case 3: The distribution of eigenfunction for a symmetric mode with 
frequency / = 110.4 kHz. ( u> = (0.0656, 0.000114) , a — 0.1 , aq = 0.033m and Re 







Fig. 4.49 Case 3: The distribution of Fourier components of streamwise velocity 
disturbance \u\j\ with Fourier modes m at different heights, (frequency / = 110.4 
kHz , w = (0.0656, 0.000114) , a = 0.1 , x x = 0.033m and Re = 531 ) 





Fig. 4.50 Case 3: The distribution of 
the azimuthal direction, (frequency / 
xj = 0.033m and Re — 531) 
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eamwise velocity perturbation profiles along 
110.4 kHz , u ; = (0.0656, 0.000114) , a = 0.1, 











b) X 3 — 0.274 mm 


Fig. 4.52 Case 4: The distribution of eigenfunction for a symmetric mode with 
frequency / = 107.6 kHz. ( uj = (0.0639, 0.00094) , a = 0.1, ii = 0.033m, and Re 
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Fig. 4.53 Case 4: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency / = 107.6 
kHz , u) = (0.0639, 0.00094) , a = 0.1 , x\ — 0.033m and Re — 531 ) 
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Fig. 4.55 Case 5: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency J 84. b 
kHz , u; = (0.05039, 0.00164) , a = 0.08 , xi = 0.033m and Re = 531 ) 
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Fig. 4.57 Case 5: The distribution of temperature perturbation profiles along the 
azimuthal direction, (frequency / = 84.8 kHz , ui = (0.05039, 0.00164) , a = 0.08. 
.Ci = 0.3505m and Re = 531) 
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a) The distribution velocity disturbance ui rea i in the azimuthal direction at 
wall-normal point X 3 = 0.123mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m. 


Fig. 4.58 Case 6 : The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 89.7 kHz. ( lo = (0.0532, 0.0023) , a = 0.08 , 
.Ti = 0.033m and Re = 531) 





Fig. 4.59 Case 6 : The distribution of disturbance profiles along the azimuthal 
direction, (frequency / = 89.7 kHz , u; = (0.0532, 0.0023) and a = 0.08, = 0.033m 

and Re = 531) 
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Fig. 4.61 Case 7: The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m at different heights, (frequency / = 80.5 
kHz , oj = (0.04785, 0.00137) , a = 0.08 , x t = 0.033m and Re = 531 ) 
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a) The distribution velocity disturbance ii\ rea t in the azimuthal direction at 
wall-normal point x z = 0.234mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity disturbance 
|u!/| with Fourier modes m at different heights x 3 . 

Fig. 4.62 Case 1: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 27.7 kHz , uj — (0.0387,0.0042 1) , cv — 0.07 , 
,Cj = 0.1978m and Re = 1350) 
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b) ^1 real 


Fig. 4.63 Case 1: Contour plot of streamwise velocity disturbance. ( frequency 
/ = 27.7 kHz , ut =(0.0589,0.00301) , a = 0.07 , x x = 0.1978m and Re = 1350 ) 





b) £3 = 0.801mm 

Fig. 4.64 Case 2: The distribution of eigenfunction for a symmetric mode with 
frequency / = 42.2 kHz. ( u> = (0.0589 ,0.00301) , a = 0.09, £1 = 0.1978m and Re 








Fig. 4.65 Case 2: The distribution of Fourier components of streamwise velocity 
disturbance with Fourier modes m. (frequency / = 42.2 kHz , u = (0.0589, 
,0.00301) , a = 0.09 , x x = 0.1978m and Re = 1350 ) 
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Fig. 4.66 Case 2: The distribution of streamwise velocity perturbation profiles along 
the azimuthal direction, (frequency / = 42.2 kHz , u> = (0.0589, 0.00301) , a = 0.09, 
X! = 0.1978m and Re = 1350) 


Fig 4 67 Case 2: The distribution of temperature perturbation profiles along the az 
imuthal direction. ( u> = (0.0589, 0.00301) and a = 0.09, X, = 0.1978m, He = 1350) 










Fig. 4.69 Case 3: The distribution of Fourier components of streamwise velocity 
disturbance \uij\ with Fourier modes m. (frequency / = 41.05 kHz , ui = (0.0572, 
0.00348) , a = 0.09 , xj = 0.1978m and Re = 1350 ) 
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wall-normal point x 3 = 0.387mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m. 


Fig. 4.70 Case 4: The distribution of eigenfunction and Fourier components lor a 
symmetric mode with frequency / = 43.36 kHz , u = (0.06056, 0.00328) , a - 0.08 
. xi = 0.1978m and Re = 1350) 
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wall-normal point x 3 = 0.387mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |«i/| with Fourier modes m. 

Fig. 4.71 Case 5: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 40.35 kHz , u = (0.0563, 0.00131) , o = 0.09 
,xi = 0.1978m and Re = 1350) 
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wall-normal point x 3 = 0.387mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |«i/| with Fourier modes m. 

Fig. 4.72 Case 6: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 41.14 kHz. (u> = (0.0574, 0.00115) , a = 0.09. 
Xi = 0.1978m. and Re = 1350) 
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Fig. 4.74 Case 6 : The distribution of disturbance profiles along the azimuthal 
direction, (frequency / = 41.14 kHz , u> = (0.0574, 0.00115) , a — 0.09, Xi — 0.197bm 
and Re = 1350) 
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wall-normal point x 3 = 0.387mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |tti/| with Fourier modes m. 

Fig. 4.75 Case 7: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 39.16 kHz. (u > = (0.0546, 0.00161) , a = 0.09, 
x 1 = 0.1978m and Re = 1350) 
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a) The distribution velocity disturbance threat in the azimuthal direction at 
wall-normal point x 3 = 0.931mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |« a /| with Fourier modes m. 

Fig. 4.77 Case 1: The distribution of eigenfunction and Fourier components for a 
symmetric mode with frequency / = 17.45 kHz , u> = (0.0385, 0.00445) , oc — 0.07, 
Xl = 0.5105m and Re = 2211) 
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Fig. 4.78 Case 1: Contour plot of streamwise velocity disturbance, (frequency 
/ = 17.45 kHz , a ; = (0.0385, 0.00445) , a = 0.07 , xj = 0.5105m and Re = 
2211 ) 
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a) The distribution velocity disturbance U\ Tea i in the azimuthal direction at 
wall-normal point = 0.931mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance \u\f\ with Fourier modes m. 


Fig. 4.79 Case 2: The distribution of eigenfunction and Fourier components for an 
anti-symmetric mode with frequency / = 19.2 kHz , 1 0 = (0.0423, 0.00351 ) , a = 0.07, 
Xl = 0.5105m and Re - 2211) 
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6 ) U-lreal 


Fig. 4.80 Case 2: Contour plot of streamwise velocity disturbance, (frequency 
/ = 19.2 kHz , w = (0.0423, 0.00351) , a = 0.07 , x x = 0.5105m and Re = 
2211 ) 
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a) The distribution velocity disturbance iiireai in the azimuthal direction at 
wall-normal point ,r 3 = 0.672mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m. 

Fig. 4.81 Case 3: The distribution of eigenfunction and Fourier components for 
an anti-symmetric mode with frequency / = 19.04 kHz. (u> = (0.0419, 0.00439) , 
a = 0.07, X\ = 0.5105m and Re = 2211) 







b) .r 3 = 1.128mm 


Fig. 4.83 Case 4: The distribution of eigenfunction for an anti-symmetric mode with 
frequency / = 23.03 kHz. ( w = (0.0507, 0.00288) , o = 0.08, Xl = 0.5105m and Re 
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Fig. 4.84 Case 4: The distribution of Fourier components of st 
disturbance \u\f\ with Fourier modes m. (frequency / = 23.03 k 
0.00288), a = 0.08, X\ = 0.5105m and Re = 2211 ) 
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Fig. 4.85 Case 4: Contour plot of streamwise velocity disturbanc 
/ = 23.03 kHz, u) = (0.0507, 0.00288), a = 0.08, x 1 = 0.5105m ai 




b) £3 = 0.85mm 


Fig. 4.86 Case 5: The distribution of eigenfunction for an anti-symmetric mode with 
frequency / = 22.04 kHz. ( w = (0.0486, 0.00427) , a = 0.08 , = 0.5105m and 

Re = 2211) 
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Fig. 4.87 Case 5: The distribution of Fourier components of streamwise velocity 
disturbance |tq/| with Fourier modes m. (frequency / = 22.04 kHz, =(0.0486, 
0.00427), a = 0.08, aq = 0.5105m and Re = 2211 ) 



Fig. 4.88 Case 5: Contour plot of streamwise velocity disturbance \u x \. ( frequency 
/ = 22.04 kHz, w = (0.0486, 0.00427), a = 0.08, x x = 0.5105m and Re = 2211 ) 
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a) The distribution velocity disturbance Ui rea / in the azimuthal direction at 
wall-normal point x 3 = 0.931mm, where it has maximum amplitude 



0 20 40 60 80 100 
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b) The distribution of Fourier components of streamwise velocity 
disturbance |ui/| with Fourier modes m. 

Fig. 4.89 Case 6: The distribution of eigenfunction and Fourier components for 
an anti-symmetric mode with frequency / = 23.45 kHz. (a; = (0.0517, 0.00411) , 
a = 0.08, xi = 0.5105m and Re = 2211) 





Fig. 4.90 Case 3: Contour plot of streamwise velocity disturbance |ui|. (frequency 
/ = 23.45 kHz , u? = (0.0517, 0.00411), a = 0.08 , x x = 0.5105m and Re = 2211 ) 




CHAPTER 5 


CONCLUSIONS 

A program is developed to investigate the linear stability of three-dimensional 
compressible boundary layer flows over bodies of revolutions. The problem is for- 
mulated as a 2-D eigenvalue problem incorporating the meanflow variations in the 
normal and azimuthal directions. Thereby, the normal mode solutions are sought in 
the whole plane perpendicular to the axial direction rather than in a line normal to 
the wall as is done in the classical theory. The case of a supersonic flow over a sharp 
cone with 5° half-included angle at 2 U angle of attack was considered. The stability 
computations were done using ID and 2D eigenvalue methods. In the case of 2D 
eigenvalue computations Implicitly Restarted Arnoldi Method was used to perform 
global eigenvalue search and those values were used as guess values for the local 2D 
eigenvalue computations. 

In the first chapter the fundamentals of the linear stability was reviewed. 
The nature of the instability in compressible and incompressible two- dimensional 
and three-dimensional boundary layers and the formulations of the stability problems 
as temporal and spatial problems were explained. Also the general historical review 
of the research on the stability were mentioned. 

Starting from the mathematical formulation of the stability problem of 
the three-dimensional compressible boundary layer in generalized curvilinear coordi- 
nate system, the numerical method and solution procedures for the case of a three- 
dimensional boundary layer over a sharp cone at an angle of attack were discussed in 
chapter2. The problem was formulated as ID and 2D eigenvalue problem. In chapter 
3, the results of meanflow computation over a sharp cone with 5° half-included angle 
at 2° angle of attack obtained using TLNS3D were presented. 
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In chapter 4 , the results of linear stability computations for the three- 
dimensional compressible boundary layer over a sharp cone at an angle of attack as 
ID and 2D eigenvalue method were presented. The stability computations were done 
at four different stations along the stream wise direction. In ID eigenvalue problem 
the stability computations were performed as in classical theory - neglecting the 
variation of the meanflow in azimuthal and streamwise direction. The ID stability 
results showed that the boundary layer region up to x\ = 0.03m is stable and most 
amplified temporal amplification rate occur around x 2 = 90°. Also the effect of cross- 
flow component was noticed to be in the middle region in the azimuthal direction and 
it manifested itself as increase in the temporal amplification rate around x 2 = 90° for 
negative mode numbers m and decrease in the temporal amplification rate for positive 
mode numbers. Unlike the ID stability results, where the eigenvalue spectrum for a 
given azimuthal mode number m consisted only a few sparse unstable eigenvalues, 
the 2D temporal eigenvalue method showed that in the 2D eigenvalue spectrum the 
unstable eigenvalues were clustered in the complex uj plane. Further, the distributions 
of eigenfunctions along the azimuthal direction were found to clustered in confined 
regions and were approximately symmetric about x 2 = 90°. Due to the huge memory 
requirement it is not possible to increase the number of points and the number of 
modes. 


5.1 Recommendations for the future work 

The major barriers in the stability computations as 2D eigenvalue method are re- 
quired memory and the CPU time. These limitations on the computational resources 
make the number of Fourier modes and the grid points in the wall normal directions 
kept to a range. Also, as explained earlier the spatial stability computations need 
much more memory resources than that of temporal eigenvalue computations. 

Thus, the recommendations for the future research are as follows; 



The stability computations of the three-dimensional compressible boundary layers 
as spatial eigenvalue problem. 

Studying of evolution of each stability mode using PSE methods. 

Computation of N-factor and prediction of transition. 

Introducing a controlled disturbance , such as a point source of the form e and 
studying the evolution of it. 
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wall-normal point x 3 = 0.931mm, where it has maximum amplitude 



b) The distribution of Fourier components of streamwise velocity 
disturbance |m/| with Fourier modes m. 


Fig. 4.89 Case 6: The distribution of eigenfunction and Fourier components tor 
an anti-symmetric mode with frequency / = 23.45 kHz. (w = (0.0517, 0.00411) , 
a = 0.08, X\ = 0.5105m and Re = 2211) 





Fig. 4.90 Case 3: Contour plot of streamwise velocity disturbance |ui|. (frequency 
/ = 23.45 kHz , u) = (0.0517, 0.00411), a = 0.08 , x\ — 0.5105m and Re = 2211 ) 
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CONCLUSIONS 


A program is developed to investigate the linear stability of three-dimensional 
compressible boundary layer flows over bodies of revolutions. The problem is for- 
mulated as a 2-D eigenvalue problem incorporating the meanflow variations in the 
normal and azimuthal directions. Thereby, the normal mode solutions are sought in 
the whole plane perpendicular to the axial direction rather than in a line normal to 
the wall as is done in the classical theory. The case of a supersonic flow over a sharp 
cone with 5° half-included angle at 2° angle of attack was considered. The stability 
computations were done using ID and 2D eigenvalue methods. In the case of 2D 
eigenvalue computations Implicitly Restarted Arnold! Method was used to perform 
global eigenvalue search and those values were used as guess values for the local 2D 

eigenvalue computations. 

In the first chapter the fundamentals of the linear stability was reviewed. 
The nature of the instability in compressible and incompressible two- dimensional 
and three-dimensional boundary layers and the formulations of the stability problems 
as temporal and spatial problems were explained. Also the general historical review 
of the research on the stability were mentioned. 

Starting from the mathematical formulation of the stability problem of 
the three-dimensional compressible boundary layer in generalized curvilinear coordi- 
nate system, the numerical method and solution procedures for the case of a three- 
dimensional boundary layer over a sharp cone at an angle of attack were discussed in 
chapter2. The problem was formulated as ID and 2D eigenvalue problem. In chapter 
3. the results of meanflow computation over a sharp cone with 5° half-included angle 
at 2 ° angle of attack obtained using TLNS3D were presented. 
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In chapter 4 , the results of linear stability computations for the three- 
dimensional compressible boundary layer over a sharp cone at an angle of attack as 
ID and 2D eigenvalue method were presented. The stability computations were done 
at four different stations along the streamwise direction. In ID eigenvalue problem 
the stability computations were performed as in classical theory - neglecting the 
variation of the meanflow in azimuthal and streamwise direction. The ID stability 
results showed that the boundary layer region up to x t = 0.03m is stable and most 
amplified temporal amplification rate occur around x 2 = 90°. Also the effect of cross- 
flow component was noticed to be in the middle region in the azimuthal direction and 
it manifested itself as increase in the temporal amplification rate around x 2 = 90° for 
negative mode numbers m and decrease in the temporal amplification rate for positive 
mode numbers. Unlike the ID stability results, where the eigenvalue spectrum for a 
given azimuthal mode number m consisted only a few sparse unstable eigenvalues, 
the 2D temporal eigenvalue method showed that in the 2D eigenvalue spectrum the 
unstable eigenvalues were clustered in the complex a; plane. Further, the distributions 
of eigenfunctions along the azimuthal direction were found to clustered in confined 
regions and were approximately symmetric about x 2 = 90°. Due to the huge memory 
requirement it is not possible to increase the number of points and the number of 

modes. 

5.1 Recommendations for the future work 

The major barriers in the stability computations as 2D eigenvalue method are re- 
quired memory and the CPU time. These limitations on the computational resources 
make the number of Fourier modes and the grid points in the wall normal directions 
kept to a range. Also, as explained earlier the spatial stability computations need 
much more memory resources than that of temporal eigenvalue computati 
Thus, the recommendations for the future research are as follows; 



The stability computations of the three-dimensional compressible boundary layers 
as spatial eigenvalue problem. 

Studying of evolution of each stability mode using PSE methods. 

Computation of N-factor and prediction of transition. 

Introducing a controlled disturbance , such as a point source of the form e ‘ m9 and 
studying the evolution of it. 
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